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PREFACE 


The  investigations  reported  herein  comprise  a  collection  of  pertinent 
literature  applicable  to  one-dimensional  stress  wave  propagation  in  soils. 
The  work  was  performed  in  connection  with  research  on  propagation  of 
ground  shock  through  earth  media  being  conducted  by  personnel  of  the 
Soils  Division,  WES,  for  DASA. 

The  work  was  accomplished  during  the  period  February  1968  through 
August  1968  under  the  supervision  of  Dr.  J.  S.  Zelasko  and  the  general 
direction  of  Messrs.  W.  J.  Turnbull,  A.  A.  Maxwell,  R.  W.  Cunny,  and 
J.  G.  Jackson,  Jr.,  of  the  Soils  Division. 

This  report  was  prepared  by  Mr.  B.  Rohani  and  the  material  contained 
herein  was  submitted  as  a  technical  paper  in  partial  fulfillment  of  the 
requirements  administrated  in  a  Civil  Engineering  graduate  course  at 
Texas  A&M  University,  College  Station,  Texas. 

COL  John  R.  Oswalt,  Jr.,  CE,  and  COL  Levi  A.  Brown,  CE,  were 
Directors  of  the  WES  during  this  investigation.  Mr.  J.  B.  Tiffany 
was  Technical  Director. 
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This  1999  report  is  a  reprint  of  a  technical  paper  prepared  by  the  author  in  1968  in  partial 
fulfillment  of  the  requirements  in  a  graduate  course  offered  by  the  Department  of  Civil  Engineering,  Texas 
A&M  University,  College  Station,  Texas.  The  work  was  performed  in  connection  with  research  on 
propagation  of  ground  shock  through  earth  media  conducted  by  members  of  the  staff  of  the  Structures 
Laboratory  (SL),  U.S.  Army  Engineer  Waterways  Experiment  Station  (WES),  Vicksburg,  MS,  a  complex 
of  five  laboratories  of  the  Engineer  Research  and  Development  Center  (ERDC),  for  the  Defense  Threat 
Reduction  Agency,  Alexandria,  VA. 

The  work  was  accomplished  by  Dr.  Behzad  Rohani,  Geomechanics  and  Explosion  Effects 
Division  (GEED),  SL,  WES,  during  February  through  August  1968  under  the  supervision  of  Dr.  J.  S. 
Zelasko  and  the  general  direction  of  Messrs.  W.  J.  Turnbull,  A.  A.  Maxwell,  R.  W.  Cunny,  and  Dr.  J.  G. 
Jackson,  Jr.,  GEED.  Dr.  Bryant  Mather  was  Director,  SL,  during  the  publication  of  this  report. 

The  purpose  of  the  study  was  to  assemble  pertinent  theoretical  literature  applicable  to  one¬ 
dimensional  stress  wave  propagation  in  soils,  rewrite  the  mathematics  in  detail  using  consistent  notations 
and  terminology,  computerize  each  solution,  prepare  instructions  for  each  computer  program,  and  to  make 
comparative  studies  with  the  various  mathematical  models  for  a  wave  propagation  problem.  The 
theoretical  developments  were  documented  in  a  format  that  can  be  used  for  engineering  training  and  self- 
study  and  are  published  for  such  purposes.  The  computer  programs,  based  on  the  technology  of  the  1960's, 
are  published  here  for  historical  and  reference  purposes. 

Commander  of  ERDC  during  the  publication  of  this  report  was  COL  Robin  R.  Cababa,  EN.  This 
report  was  published  at  the  WES  complex  of  ERDC. 
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Introduction 


A.  Background 


The  first  major  research  concerning  the  behavior  of  soils  under 
transient  loadings  was  sponsored  by  the  United  States  Army  Corps  of 
Engineers  to  aid  in  the  design  of  underground  protective  structures.  The 
work  was  carried  out  at  the  Massachusetts  Institute  of  Technology  (MIT) 
under  the  directions  of  D.  ¥.  Taylor  and  R.  V.  Whitman.  A  brief  summary 
of  this  work  was  published  by  Whitman  (1)  in  1957;  it  states 

A  hydraulic  apparatus  and  special  instrumentation 
were  constructed  to  test  triaxial  soil  samples.  Failure 
was  achieved  in  times  as  short  as  0.001  second.  Curves 
of  compressive  strength  versus  strain  rate  (rapidity  of 
loading)  were  determined  for  cohesive  soils,  dry  sand, 
and  saturated  sands.  Transient  pore-water  pressures  were 
recorded  during  tests  on  saturated  sands. 

Another  apparatus  was  constructed  to  study  wave 
propagation.  Soil  samples,  2  in.  in  diameter  and  32  in. 
long,  were  struck  at  one  end  by  a  ram.  Results  for  a 
dry  sand  were  compared  with  theoretical  solutions  for 
the  wave  propagation  problem. 

Other  tests  were  devised  to  study  creep  and  re¬ 
laxation  phenomena  in  dry  sand,  and  to  study  the  permea¬ 
bility  of  saturated  sands  to  pressure  gradients  applied 
suddenly. 

Whitman  compared  the  results  of  the  wave  propagation  experiments 
with  the  one-dimensional  rate -independent  plastic  wave  propagation  theory 
developed  by  Von  Karman  and  Duwez  (see  Chapter  II,  Section  ~D).  The 
theory  did  not  predict  the  initial  peak  stress  or  spike  at  the  impact  end 
of  the  sand  column  (fig.  1-1)  and  he  concluded  that  lateral  inertia  effects 
were  responsible  for  the  occurrence  of  the  stress  peak.  In  1961,  B.  R. 
Parkin  (2)  developed  a  rate-dependent  elastic-plastic  theory  to  study 
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one-dimensional  stress  wave  propagation  in  sand.  This  theory  predicted 
the  occurrence  of  the  impact  stress  peak  of  the  MIT  experiments  without 
considering  lateral  inertia  effects  (fig.  1-1).  Parkin  concluded  that 
the  initial  peak  stress  was  due  entirely  to  interaction  of  strain-rate 
sensitivity  and  static  stress-strain  properties  of  the  medium  with  the 
elastic  compliance  of  the  impacting  stress  gage.  Parkin's  conclusion 
aroused  much  interest  and  resulted  in  a  symposium  (3)  on  impact  waves 
in  sand. 

The  subject  of  wave  propagation  in  soils  was  given  much  attention 
after  Parkin's  work.  New  advances  in  nuclear  technology  and  space  ven¬ 
tures,  in  particular,  provided  considerable  Impetus  towards  this  effort. 
Presently  (1968)  a  number  of  research  organizations  are  working  on  dif¬ 
ferent  phases  of  the  problem  with  an  overall  objective  of  predicting  wave 
propagation  phenomena  in  an  in  situ  soil  mass. 

B.  The  problem 

The  problem  to  be  considered  in  wave  propagation  in  soils  is  that  of 
predicting  the  form  and  the  effect  of  an  input  wave  after  it  has  propagated 
through  the  soil,  at  a  specified  point  in  space  and  time.  In  particular, 
the  attenuation  of  peak  stress  and  particle  velocity  with  depth  is  of 
interest.  The  solution  of  this  problem,  like  any  other  boundary  value 
problem,  requires  a  knowledge  of: 

a.  the  equations  governing  the  motion  of  the  soil  mass, 

b.  the  equation  of  continuity  which  expresses  the  conservation  of 
mass, 

c.  the  constitutive  equations  which  relate  stress  to  strain  and 
strain  rate,  and 

d.  the  boundary  and  initial  conditions. 
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The  solution  of  the  equations  of  motion  with  specified  boundary  and  initial 
conditions  and  constitutive  equations  is  the  solution  to  the  problem. 

Major  difficulties  arise  when  attempts  are  made  to  solve  such  problems  for 
soils.  The  fact  that  the  stress- strain  relations  for  soils  are  highly 
non-linear,  hysteretic  and  rate-dependent  makes  the  task  of  obtaining  an 
analytical  solution  a  difficult  one.  The  problem  becomes  even  more 
complicated  when  one  considers  the  inhomogeneity  of  in  situ  soil  masses 
and  the  fact  that  a  soil  medium  is  not  strictly  a  continuum. 

C.  Equations  of  motions  and  continuity 

The  equations  of  motion  expressing  the'  equality  of  applied  force  and 
time-rate  of  change  of  momentum,  neglecting  body  forces,  are  given  as 


do 


xx 


dx 


do 


xy 


By 


do 

xz 

dz 


P 


dt 


do 


yx 


do.. 


dx 


yy 


do 


By 


yz 


dz 


1-1 


P 


do 


zx 


dx 


do 


zy 


By 


do 


zz 

dz 


where 

p  is  the  mass  density  of  the  material, 
ux  ,  Uy  ,  uz  are  the  velocity  components  in  the  x  , 
y  ,  and  z  Cartesian  directions 
respectively,  and 

ct  ,  o  ,  o  ,  o  ,  etc.,  are  components  of  the  stress  tensor  in 
xx  5  xy  5  xz  5  yx 

Cartesian  coordinates. 

These  equations  of  motion  will  hold,  irregardless  of  the  stress-strain 
behavior  of  the  medium. 
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The  equation  of  continuity  expressing  the  conservation  of  mass  in 
Cartesian  coordinates  is  given  by 


d£ 

dt 


+ 


1-2 


Equations  1-1  and  1-2  are  4  equations  in  10  unknowns  (6  stress  components, 

3  velocity  components,  and  the  mass  density  of  the  medium).  To  solve 
Equations  1-1  and  1-2  six  additional  equations  are  needed.  These  addi¬ 
tional  equations  must  be  obtained  from  the  constitutive  relations  of  the 
material . 

D.  Constitutive  relations 

Constitutive  relations  are  equations  which  relate  in  some  manner  the 
respective  components  of  the  stress  and  strain  tensors.  These  relations 
must  be  determined  from  property  tests  conducted  on  samples  of  the  material 
of  interest  to  measure  the  appropriate  material  properties.  For  a  lin¬ 
early  elastic,  isotropic  and  homogeneous  material  only  two  constants  are 
needed  to  completely  formulate  the  constitutive  equations.  These  constants 
can  be  determined  by  performing  appropriate  property  tests.  The  most 
conventional  property  test  is  the  simple  tension  test  from  which  both 
Young's  modulus  of  elasticity,  E,  and  Poisson's  ratio,  v,  can  be  obtained. 
These  constants  are  related  to  other  elastic  material  parameters  as 
follows : 


M  = 


E(1  -  v)  3K(1  -  v)  2G(1  -  v) 

(1  +  v)(l  -  2v)  (1  +  v)  (l  -  2v) 


1-3 


where 

M  is  the  constrained  modulus  of  elasticity 
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K  is  the  hulk  modulus  of  elasticity 

G  is  the  shear  modulus 

The  constitutive  equations  for  a  linearly  elastic,  isotropic  and 
homogeneous  material  can  he  written  in  terms  of  any  two  of  these  elastic 
constants.  In  terms  of  Young' s  modulus  and  Poisson's  ratio,  these 
equations  are  given  hy 

°xx  ~~  (1  -  2v)(l  +  v)  exx  +  v  ^eyy  +  ezzO 

°yy  =  (1  -  2v) (1  +  v)  eyy  +  V  ^exx  +  ®zzi] 

CTzz  ~  (1  -  2v)  (1  +  v)  ~  ezz  +  v  ^exx  +  eyy^-l 

1-4 

_  E 

°xy  2(l  +  v)  xy 

_  E 

°xz  2(1  +  v)  sxz 

E 

yz  2(1  +  v)  yz 

where 

e  ,  e  ,  e  ,  e  ,  etc.,  are  components  of  the  strain  tensor  in 

xx  5  xy  5  xz  3  yx  ’  J 

Cartesian  coordinates. 

Stress-strain  relations  for  soils  are  not  unique  and  in  general  they 
assume  a  variety  of  forms  depending  upon  many  factors  such  as  state  of 
stress,  previous  stress  history,  rate  of  loading,  degree  of  satura¬ 
tion,  etc.  Unlike  the  linear  elastic  material,  each  property  test  for 
soils  will  yield  a  highly  non-linear  relation  between  the  appropriate  stress 
and  strain  components.  There  are  no  representative  equations  that  can 
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completely  relate  these  components  in  a  rigorous  manner  as  is  done  in  the 
classical  theory  of  linear  elasticity  (equation  1-4). 

In  order  to  minimize  rate  effects,  soil  property  tests  for  dynamic 
problems  should  be  conducted  at  rates  appropriate  to  the  type  of  dynamic 
problem  investigated.  Two  types  of  dynamic  soil  property  tests  currently 
in  use  are  the  uniaxial  strain  test  (fig.  1-2)  and  the  triaxial  compression 
test  (fig.  1-3).  Radial  symmetry  is  established  in  both  of  these  tests. 

The  slope  of  the  uniaxial  strain  test  defines  the  constrained  modulus  M 
while  the  slope  of  the  triaxial  compression  test  curve  defines  the  shear 
modulus  G  where  radial  strain  is  known.  The  complexity  of  formulating 
constitutive  equations  for  soils  is  illustrated  by  the  diverse  nature 
property  test  results. 

E.  Boundary  loads 

The  dynamic  boundary  loads  which  are  of  interest  in  wave  propagation 
studies  in  soils  are  shown  in  fig.  1-4.  These  loads  may  be  generated  by 
explosions  or  vibrating  machinery.  Fig.  I-4a  depicts  a  typical 
overpressure-distance  curve  resulting  from  a  nuclear  explosion.  This  type 
of  loading  is  probably  the  most  complicated  boundary  load  that  has  to  be 
dealt  with  in  solving  a  dynamic  problem.  The  fact  that  the  velocity  of 
the  traveling  shock  wave  changes  as  it  moves  away  from  ground  zero  is  a 
major  component  of  the  complexity  of  this  problem. 

Fig.  I-4b  shows  a  single  pulse  characterized  by  a  peak  stress  at  the 
shock  front  and  an  exponentially  decaying  behavior  thereafter.  This  type 
of  loading  may  be  generated  by  a  single  charge  explosion. 

A  steady- state  type  input  generated  by  vibrating  machinery  is  shown 
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Fig.  1-2.  Typical  load-unload  stress-strain  curve  for  soils  in  uniaxial  strain. 
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Fig.  1-4.  Typical  boundary  loads. 
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in  fig.  I-4c.  This  type  of  loading  is  usually  of  low  amplitude  and  rela¬ 
tively  long  duration. 

F.  Scope  of  the  present  study 

In  this  study,  the  response  of  a  semi-infinite  homogeneous  body  of 
soil  due  to  a  time -dependent  pressure  wave  applied  at  the  free  surface  of 
the  medium  is  to  be  studied.  The  loading  is  assumed  to  be  applied  uni¬ 
formly  over  the  entire  surface  of  the  medium.  Under  this  assumption,  one- 
dimensional  geometry  is  obtained,  i.e.,  there  sire  no  lateral  deformations 
and  the  equations  of  motion  (equation  1-1)  and  continuity  (equation  1-2) 
simplify  as  follows; 


Su  _  So 
"  dt  -  dz 


1-5 


|f  +  pp.  0 

dt  dz 


1-6 


The  subscripts  have  been  dropped  since  only  vertical  components  are 
to  be  considered. 

The  stress-strain  relation  required  by  the  one-dimensional  geometry 
is  the  uniaxial  strain  relation  (fig.  1-2).  This  relation  is  not  unique 
and  there  are  a  number  of  regions  in  this  curve  in  which  different  types 
of  material  behavior  are  dominate.  For  example,  for  dry  cohesionless 
materials,  during  the  early  stages  of  loading,  the  stress-strain  curve  may 
be  concave  downward  reflecting  a  rearrangement  of  particles.  As  the  load 
is  increased  the  particles  begin  to  lock  together  and  the  strains  that 
take  place  are  due  primarily  to  the  deformation  of  the  particles  at  the 
points  of  contact.  As  the  stress  continues  to  increase,  the  contact  forces 
become  so  large  that  the  particles  begin  to  crush.  As  the  stress  is 


11 


further  increased,  the  particles  lock  again  and  the  curve  becomes  concave 
to  the  stress  axis  with  the  particles  smaller  and  more  angular  than 
before  re  stilting  in  a  progressively  higher  modulus.  During  unloading, 
the  curve  exhibits  a  higher  modulus  than  the  loading  curve  and  permanent 
strain  occurs.  For  cohesive  soils,  rate -dependency  is  also  a  factor  to 
be  considered  in  utilizing  the  stress- strain  relation  of  fig.  1-2. 

To  overcome  the  mathematical  difficulties  in  the  solution  to  these 
problems,  real  soil  stress- strain  response  may  be  idealized  by  various 
linear  and/or  nonlinear  hysteretic  approximations  as  shown  in  fig.  1-5* 

A  considerable  body  of  scientific  literature  on  one- dimensional  stress 
wave  propagation  for  such  models  has  been  published  in  recent  years  by 
various  researchers  both  in  the  United  States  and  abroad.  It  was  the  aim 
of  this  study  to  assemble  a.1 1.  pertinent  literature  applicable  to  wave 
propagation  in  soils,  rewrite  the  mathematics  in  detail,  computerize  each 
solution,  prepare  operating  instructions  for  each  computer  program,  and  to 
make  comparative  studies  with  the  various  models  for  a  wave  propagation 
problem.  In  Chapter  II,  the  elements  of  elastic,  plastic,  visco-elastic 
and  shock-wave  propagation  are  reviewed  in  detail  as  a  background  for  the 
remainder  of  the  report. 
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Nonlinear  hysteretic 


Fig.  1-5.  Idealized  stress-strain  curves  for  soils 
in  uniaxial  state  of  strain. 
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II 


Elements  of  One -Dimensional 
Stress  Wave  Propagation 


A.  Historical  background 

Interest  in  the  theory  of  propagation  of  longitudinal  waves  dates 
frcm  the  17th  century  when  Newton  attempted  to  solve  the  problem  of  sound 
wave  propagations,  a  solution  finally  completed  by  Lord  Rayleigh  in  the 
19th  century  (4 ) .  Much  of  the  theory  in  this  period  was  developed  for 
materials  with  highly  idealized  stress-strain  relations.  Hie  theory  of 
propagation  of  longitudinal  waves  in  elastic  bars  is  presented  in  a 
number  of  publications  (5 ) .  The  solution  for  this  case  is  virtually 
complete  if  the  lateral  inertial  forces  are  negligible.  ‘  Pochhammer  solved 
the  problem  considering  lateral  inertia  for  a  cylindrical  bar  (6).  Stress 
wave  propagation  in  visco-elastic  materials  has  been  studied  by  Kolsky  (5) 
and  Morrison  (7)  for  linear  visco-elastic  materials.  An  excellent  summary 
of  one-dimensional  wave  propagation  theory  for  standard  models  of  visco¬ 
elasticity  as  well  as  three  and  four  parameters  models  is  given  by  Kolsky 
(1963).  The  first  development  of  plastic  wave  propagation  theory  is 
credited  to  Donnell  (8).  He  suggested  that  if  the  stress-strain  curve  for 
a  plastic  material  could  be  approximated  by  two  straight  lines,  then  the 
time  history  of  a  wave  propagating  through  that  material  can  be  deduced  by 
superposition  of  the  time  histories  of  the  two  elastic  waves.  The  com¬ 
plete  theory  of  one-dimensional  plastic  wave  propagation  was  developed 
during  the  Second  World  War.  Solutions  were  obtained  independently  by 
Von  Karman  (9),  Taylor  (10) ,  and  Rokhmatulin  (ll).  The  theory  was  later 
extended  by  White  and  Griffis  (12)  to  handle  both  plastic  and  shock  wave 
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propagation  and  by  Malvern  to  handle  materials  which  exhibit  strain  rate 
effects  (13) . 

B.  Longitudinal  elastic  waves 

In  the  simple  theory  of  longitudinal  wave  propagation  (no  lateral 
strain)  it  is  assumed  that  plane  cross  sections  remain  plane;  only  axial 
stresses  are  considered,  being  uniformly  distributed  over  the  cross  section 
(fig.  II-l) .  The  cross-sectional  area  of  the  element  is  denoted  by  A  , 
the  thickness  by  dz  and  the  mass  density  by  p  .  The  displacement  of  the 
element  in  z  direction  is  given  by  u  .  Newton's  second  law  of  motion 
for  the  element  of  fig.  II-l  gives 


-p  A 


II-l 


or 


u  d 
St2 


dz 


II-2 


where  a  is  the  axial  compressive  stress.  For  infinitesimal  strain,  the 
equation  of  continuity  in  the  one-dimensional  case  reduces  to  the  defini¬ 
tion  of  strain 


II-3 


where  e  is  the  axial  compressive  strain.  The  constitutive  relation  for 
an  elastic  medium  with  no  lateral  motion  allowed  is  given  by 


a  =  M  e 
o 


II-4 


where  M  is  the  constrained  modulus  of  elasticity.  M  is  related  to 
o  o 
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P(t) 


Fig.  II-l.  One-dimensional  wave  propagation. 
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Lame's  constants  X  and  p  or  to  Young's  modulus  E  and  Poisson's 
ratio  v  as  follows 


M  =  X  +  2p 
o 


E(1  -  V) 

(1  +  v)(l  -  2v) 


II-5 


Substituting  equation  II-3  in  equation  II -4  and  differentiating  with  re¬ 
spect  to  the  space  variable  results  in 


8z 


II -6 


The  substitution  of  equation  II-6  into  equation  II-2  yields  the  well 
known  equation 


d2u  p2  d2u 

3t2  =  0  Bz2 


II -7 


where  Cq  is  the  propagation  velocity  of  longitudinal  waves 


p(l  +  V)(l  -  2v) 


II-8 


The  general  solution  of  equation  II -7  is  of  the  form 


II -9 


This  solution  represents  the  sum  of  two  traveling  waves.  The  first  term 
represents  a  wave  traveling  in  the  positive  z  direction,  and  the  second 
term  represents  a  wave  traveling  in  the  negative  z  direction.  For  a 
medium  of  infinite  extent  where  the  waves  travel  only  in  the  positive 
z  direction  the  solution  is  given  by 


11-10 
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Differentiating  both  sides  of  equation  11-10  with  respect  to  z  and  t  , 
respectively,  gives 


Su 

dt 


11-11 


11-12 


where  f'  denotes  differentiation  with  respect  to  the  argument  (-*)• 
is  the  particle  velocity  u  .  Combining  equations  11-11  and  11-12  gives 


du  _  _  Su 

£t  o  Sz- 


Using  the  definition  of  strain  (equation  II-3) 


u  =  C  e 
o 


Substituting  for  strain  from  equation  II -4 


M 

O  *  r~i  • 

0  =  —  u  = 

O 


11-13 


II-14 


11-15 


Equation  U-l4  relates  particle  velocity  and  strain;  equation  11-15  re¬ 
lates  particle  velocity  and  stress. 

The  function  <*  -  ft)  must  be  evaluated  from  the  boundary  con¬ 


ditions.  At  z  =  0 


a(0,t)  =  P(t) 


11-16 


where  P(t)  is  a  known  arbitrary  input  pressure.  A  convenient  time  param¬ 
eter,  defined  by 


5 


is  usually  used  to  adjust  the  real  time  at  the  surface,  by  the  factor  — 

co 

to  account  for  the  travel  time  of  the  wave  to  the  point  of  interest. 
Substituting  equation  11-12  for  particle  velocity  into  equation  II -15 
and  utilizing  equation  11-17  gives 


0  =  <V'(t  -§;)= 

conditions  at  the  boundary,  equation  II -16,  require  that 

PCQf'(t*)  =  P(t*) 


11-18 


11-19 


where  P(t*)  is  the  surface  overpressure  at  time  t  =  t*  .  An  expression 


for 


K*  ■  fe) 


can  then  be  derived  from  equation  11-19 . 

The  expressions  for  stress,  strain,  and  particle  velocity,  in  terms 
of  the  surface  overpressure  are  given  by 


a  =  P(t*) 

e  =  P(t*)/pC^  11-20 

and 

u  =  P(t*)/pCQ 


Equation  11-20  demonstrates  that  there  is  a  linear  relation  between  the 
three  given  wave-form  parameters  and  the  surface  overpressure  P(t)  , 
and  that  for  conditions  of  one-dimensional  wave  propagation  in  an  elastic 
homogeneous  medium  an  unaltered  wave  form  propagates'  through  the 
medium. 

The  acceleration  at  a  point  -can  be  determined  from  equation  11-20  by 
differentiating  the  expression  for  particle  velocity 
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11-21 


d2u  Id 
St2  '  eco  dt 


Equation  11-21  demonstrates  that  the  acceleration  depends  on  the  time  rate 

of  change  of  the  surface  overpressure. 

The  absolute  displacement  of  a  point,  at  time,  t  ,  can  he  obtained 

by  integration  of  equation  11-20  for  particle  velocity 

t 


u 


P(t*)  dt 


11-22 


The  integral 


in  equation  11-22  is  the  input  impulse  (area  under  the  surface  overpressure- 

2 

time  curve)  between  the  times  t  =  tt-  and  t  =  t  . 

o 

C.  Longitudinal  elastic  waves  in  layered  elastic  media 

Consider  the  semi -infinite  body  of  fig.  II -2  made  up  of  two  elastic 
media  with  different  properties  separated  by  a  plane  interface.  The  prop¬ 
erties  in  the  two  layers  are  subscripted  o  and  1  .  A  wave  in  the  first 
layer  incident  to  the  interface  is  given  by  equation  11-10 


u 


11-10 


The  corresponding  incident  stress,  ,  is  given  by  equation  II-18 


c,  =  P  C  J 


o  O' 


■(•  - 1) 


11-18 


At  the  interface  a  reflected  wave  will  form  and  travel  in  the  negative  z 
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direction  in  layer  1.  The  reflected  -wave  is  given  by  the  second  term  of 
equation  II-9 


u=g(rt) 

The  reflected  stress  wave  is  given  as 

‘r  -  -  +  f) 

A  transmitted  wave 

u  =  h(t'il) 


II -9 


11-23 


11-24 


will  also  form  at  the  interface  and  will  travel  in  the  positive  z  direc¬ 
tion  in  layer  2.  The  transmitted  stress  wave  is  given  as 


°t  PlClh 


■(*■« 


11-25 


The  corresponding  particle  velocities  for  the  incident,  the  reflected  and 
transmitted  waves,  from  equation  II -18,  are  respectively 

o. 


'f-rVirf 
\  0/0  0 

e'(t  +  §-V$- 
\  0/  0  0 


II -26 


11-27 


■  y  •  °a 


11-28 


There  are  two  distinct  conditions  which  must  be  satisfied  on  both 
sides  of  the  interface 
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(a)  For  equilibrium 


c. 

x 


+  a 


r 


CJ 


t 


(b)  For  compatibility  of  particle  velocity 


P  C 
o  o 


P  C 
o  o 


Plcl 


11-29 


11-30 


The  negative  sign  in  front  of  the  reflected  particle  velocity  is  due  to 
the  fact  that  its  direction  is  in  the  negative  z  direction.  Substituting 
for  cr  =o,-c.  from  equation  11-29  in  equation  11-30  and  solving  for 

I*  Xf  X 

0.  gives 


o,  =  2o. 
t  1 


P1C1 


PC  +  PC 
00  11 


11-31 


Substituting  for  o t  in  equation  11-31  from  equation  11-29  and  solving 
for  or  gives 


o 


r 


P1C1  -  poCo\ 

PcCo  +  plCl  j 


11-32 


Interesting  conclusions  can  be  drawn  from  equations  11-31  and  11-32  re¬ 
garding  the  nature  of  the  reflected  and  transmitted  waves  at  the  inter¬ 
face.  For  instance,  at  a  free  surface,  p^C^  =  0  .  From  equation  11-31, 
o^  =  0  ;  there  is  no  transmitted  wave.  From  equation  11-32,  o^  =  -ck  ; 
the  reflected  wave  is  equal  in  magnitude  to  the  incident  wave,  but  opposite 
in  sign.  Thus,  a  compressive  wave  is  reflected  as  a  tensile  wave  at  a  free 
boundary.  At  a  rigid  boundary,  p^C^  =  00  .  From  equation  11-32,  the  re¬ 
flected  wave  0r  =  5  the  reflected  wave  is  equal  in  sign  and  magnitude 
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to  the  incident  wave.  Therefore,  at  a  rigid  boundary  the  stresses  are 
doubled. 

In  general,  from  equation  11-32,  if  PoCo  >  ,  a  compressive  wave 

is  reflected  as  a  tensile  wave.  If,  on  the  other  hand,  poCo  <  P1C1  ,  a 
compressive  wave  is  reflected  as  a  compressive  wave. 

D.  Propagation  of  plastic  waves  in  strain-rate  independent  media 

The  first  development  of  plastic  wave  propagation  was  published  by 
Donnell  (8)  in  1930.  Donnell  suggested  the  use  of  the  principle  of 
superposition  if  the  stress-strain  curve  for  a  plastic  material  can  be 
approximated  by  two  straight  lines.  Consider  the  plastic  stress-strain 
curve  and  its  approximation  in  fig.  II-3&.  If  a  step  pulse  of  stress  aQ 
is  applied  on  the  surface  of  a  semi-infinite  body  having  such  a  linearized 
stress-strain  curve,  the  step  pulse  will  propagate  into  the  body  with  a 
velocity 

Co  =  ^/iTTp  11-33 

as  shown  in  fig.  II-3b.  If  a  second  step  pulse  of  stress  is  super¬ 

imposed  on  the  body  it  will  propagate  with  a  velocity 

C1  =  VMi/p  11-34 

Thus,  at  some  time,  t  ,  the  front  of  the  wave  a  <  oq  would  have  traveled 
to  a  position 

=  C  t  11-35 

l  o 

while  the  remainder  of  the  wave  a  >  o  would  have  traveled  only  to  a 

o 

position 
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Stress 


Fig.  II-3.  Donnell's  approximation  of  plastic  wave  propagation 


as  shown  in  fig.  II -3b.  Donnell's  development  simply  indicates  that  if  the 
two  stress  pulses  were  applied  simultaneously,  the  action  would  have  been 
the  same  and  the  wave  front  would  have  the  appearance  shown  in  fig.  II -3b. 
Thus,  from  Donnell's  deductions,  one  would  expect  the  higher  stress  levels 
(for  the  real  plastic  stress-strain  curve)  to  travel  at  lower  velocities. 
The  gap  (C  t  -  C^t)  between  the  two  wave  fronts,  fig.  II-3b,  would  increase 
with  time  and  the  sharp-fronted  stress  pulse  would  rapidly  become  softened. 
Obviously,  stresses  greater  than  the  yield  stress  would  not  propagate  in  a 
plastic  medium.  From  the  above  reasoning,  one  would  expect  that  the  wave 
front  corresponding  to  the  simultaneous  action  of  two-step  waves  aQ  and 
o  on  the  surface  of  a  medium  having  a  stress-strain  curve  such  as  the 
one  shown  in  fig.  II-4a  would  have  the  form  shown  in  fig.  II-4b.  The 
mathematical  solution  of  plastic  wave  propagation  obtained  by  Von  Karman 
and  others  during  the  Second  World  War  verifies  the  correctness  of 
Donnell's  superposition  approach  as  will  be  shown  subsequently. 

Consider  a  semi-infinite  body  extending  from  z  =  0  to  z  =  oo  , 
having  a  unique  stress-strain  curve  which  is  plastic  in  character  or 
concave  to  the  strain  axis,  as  shown  in  fig.  II-5.  The  stress-strain 
curve  can  be  given  as 

aP  =  aP^ep)  11-37 

where  the  subscript  P  indicates  that  the  curve  is  concave  to  the  strain 
axis,  or  "plastic."  The  problem  being  one-dimensional,  the  equations  of 
motion  and  continuity  from  the  previous  section  are  given  as 


II -2 


and 


€P 


8z 


Su 

5z 


II-3 


For  any  point  on  the  stress-strain  curve,  the  constrained  modulus,  Mp  ,  is 


S°P 

MP(£P}  =  ^ 


11-38 


The  initial  tangent  modulus  is  Mq  .  From  equations  11-37,  II-3,  and  II-38 

80  80  8e  2 

sr "  ss;  § r =  -  W  gr 

Substituting  equation  11-39  into  equation  II -2  the  equation  of  motion 
becomes 

,2  ^2 


d  u  „  r  v 

at2  "  p  p  az2 


II-4o 


For  boundary  conditions,  assume  that  a  constant  velocity  -  u^  ,  corre¬ 
sponding  to  a  compressive  stress  -  0^  ,  is  suddenly  imposed  and  maintained 
at  the  surface  of  the  body  ( z  =  0 ) . 

Then 

u  =  -  Ujt  at  z  =  0 


and 


u  =  0  at  z  =  co 


To  integrate  equation  II -40,  assume  that  the  compressive  strain  e  is 

P 

a  function  of  z/t  =  f  ,  i.e. 


ep  =  4W 


II-41 
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Then  dz 
mined  as 


t  dt  and  dt  - - dt  ,  the  displacement  u  is  then 

z 


z, 

-  / 

00 


du 

57 


dz  = 


z  z. 

/ 

00  00 


fit)  dz 


t 

=  -  t  J  fit)  di|r 

oo 


d2 


u 


The  second  space  derivative,  — ^  ,  becomes 

dz 


2ji.  -€=-§-  Kt)  =  -  hM 

^2  Sz  P  ZzrKf'  t 


and 


du 

5t 


t  t  ' 

J  fit)  ^t-tg  J  fit)  ^ 

oo  oo 


But 


^  dt  ^  dt  dt  t  Sijr 


d  dt 


and 


£  d 
t  5t 


t 

J  fit)  dt  =  -  |  fit)  =  ~ 

oo 


tfit) 


Substituting  equation  11-46  into  equation  11-44 


du 

5t 


t 

=  -J 

oo 


fit)  dt  +  tf  it) 


Therefore 


deter- 


11-42 


11-43 


11-44 


11-45 


11-46 


11-47 
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d2u  _ 
St2  = 


8 

St 

r 

t  1 

+/(*)  -  /  fit)  dy 

L 

00  J 

>  .  's 

r  t 

1 

ay  o 

5t  5y 

tfit)  -  J  fit)  dy 

00 


=  -  %  [/(♦)  +  -  /(♦)] 
=  -%♦/'(♦)  =  -  T-  f'M 

t  t 


11-48 


Substituting  equations  11-43  and  11-48  into  equation  II-40  the  equation  of 
motion  becomes 


p  *-#’(*>  =Mp(€p)idll 


11-49 


or 


f'W 


Pt  -  Mp(ep) 


=  0 


11-50 


Equation  11-50  indicates  that  either 


=0  or  ^(y)  =  constant 


11-51 


or 

Mp(ep)  =  py2  H-52 

The  first  solution  given  by  equation  11-51  indicates  that  there  is  a  strain 
ep  which  is  constant.  The  displacement,  from  equation  11-42  and  the  boun¬ 
dary  conditions,  can  be  written  as 
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u  =  -  V  “  £iz 


11-53 


The  second  solution  given  by  equation  11-52  corresponds  to 


'Mp(ep) 


11-54 


The  complete  solution  of  the  problem  can  be  written  as  follows 

(a)  The  strain  ep  =  ei  =  constant  ,  from  z  =  0  to  z  =  Cpt 

where  ep  is  the  plastic  increment  corresponding  to  stress  cp  induced  by 

the  applied  velocity  u_  .  - - 

z  /V€p) 

(b)  Between  z  =  C,t  and  z  =  CQt  ,  the  relation  ^  =/J  — - - 

holds.  The  strain  varies  and  each  strain  increment  propagates  with  a 

/V£P)  FI 

velocity  A/  -  p  ■-■■■■  .  Cq  is  the  elastic  velocity  —  . 

(c)  e_  =  0  for  z  >  C  t  . 

This  strain  distribution  is  shown  in  fig.  II-6.  The  elastic  strain  front 
travels  with  velocity  Cq  and  an  amplitude  .  The  peak  plastic  strain 
front  propagates  with  a  velocity  Cp  and  a  constant  amplitude  . 

Between  the  two  fronts  the  propagation  velocity  of  any  strain  increment 
depends  upon  the  local  slope  Mp(€p)  and  the  strain  varies  from  eo  at 
the  elastic  front  to  ep  at  the  peak  plastic  front.  The  propagation 
velocity  of  the  plastic  wave  peak,  ,  can  be  related  to  ep  and  up  in 
the  following  way.  From  the  boundary  conditions  and  equation  11-42 


z.=u 

=  -  f  Ht) 


11-55 
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Strain 


Fig.  II-6.  Strain  distribution  for  a  plastic  pulse  in  a 
bar  according  to  the  strain-rate  independent 
theory . 
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ft|N 


or 


11-56 


This  integral  is  the  area  under  the  curve  in  fig.  II-6.  It  may  be  written 
as 


11-57 


Substituting  for 


f-  from  equation  11-5^ 

U 


*1 


iW 


de 


11-58 


The  relation  Mp(ep)  is  known  from  the  stress -strain  curve  for  the  medium. 
Equation  II-58  gives  the  relation  between  e  and  u^  .  The  velocity  Cp 

/Vs) 

is  equal  to  the  value  of  W  p  at  ^  .  The  stress  distribution  is 
deduced  from  the  distribution  of  strain  and  the  stress-strain  curve. 

E.  Propagation  of  stress  waves  in  linear  visco-elastic  media 

Linear  visco-elastic  materials  are  generally  defined  by  linear  dif¬ 
ferential  equations  which  relate  stresses  to  strain  and  strain  rate.  The 
elasticity  of  the  material  is  presented  by  linear  springs  while  the 
viscosity  of  the  material  is  presented  by  viscous  elements  called  dashpots. 
The  elastic  and  viscous  elements  are  combined  to  form  the  constitutive 
equations  of  the  material.  The  governing  differential  equations  describing 
the  propagation  of  one -dimensional  stress  waves  in  visco-elastic  materials 
are  obtained  from  the  equilibrium  equation  of  motion,  equation  II-2, 

3^ 


by  substituting  the  appropriate  constitutive  equation  for  stress  a 


S^u  do 


II-59 


Various  models  of  visco-elastic  materials  can  be  constructed  by  suitable 
combinations  of  springs  and  dashpots.  The  best  known  models  are  the  two- 
element  Kelvin-Voigt  and  Maxwell  solids  and  the  three-element  standard 
linear  solid  shown  in  fig.  II-7-  The  constitutive  relation  for  the  Kelvin- 
Voigt  model  is,  from  fig.  11-7} 


For  the  Maxwell  model 


a  =  Me  + 


T|  ££ 

dt 


(Kelvin-Voigt ) 


or 


a 


=  -n 


d(£  -  61> 

dt 


but 


1  do 
M  dt 


11-60 


II-61 


11-62 


11-63 


Substituting  equation  11-63  in  equation  11-62  results  in 


-  de  ~n  dg 
0  dt  ”-M  dt 


For  the  standard  linear  solid 

0  -  “isi  +  ”  3T  -  Vs 


Since  e  =  an(3- 


(Maxwell)  11-64 


11-65 
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a  a 


Fig.  II-7.  Models  for  linear  visco-elastic  materials 
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11-66 


del  de  3^_  dp 
dt  it  "  M  dt 

O 


equation  11-65  becomes 


_  de  "H  da 

a  “  ^  dt  “  M  dt  “l  1 
o 


11-67 


but 


therefore 


€1  ”  €  -  €2  £  "  M_ 


o  = 


Mo^l  _  .  de  da 


e  + 


Mq  +  M-|_  Mq  +  dt  Mq  +  ^  dt 


77-  (Standard  Linear  Solid)  11-68 


Substitution  of  equations  II-60,  11-64,  and  11-68  for  a  ,  in  equation 

11-59 j  will  result  in  the  following  equations  of  motion  for  the  Kelvin - 

Voigt,  Maxwell,  and  standard  linear  solids,  respectively  (note  that  strain 

8u 

and  displacement  are  related  by  equation  II -3,  that  is,  e  =  -  ^  ): 

d\ 


*2 

„  o  u  «  °  u  ,  T1 

p  — =  M  — -r  +  ^  _ 

dt  dz  dz  dt 


( Kelvin-Voigt )  II-69 


d2u  _ 
p  — £  =  *1 


dt£ 


d\ 


dz2dt 


p  d\ 
dt' 


M  ^3 


(Maxwell)  11-70 


2 

\  o  u  d  u  _  CJ-u 

p(Mo  +  V  V2  “  “A  V2  -  p11  V3  + 


A  .  S3u 


dt 


dz£ 


dt' 


dz  dt 


(Standard  Linear  Solid)  11-71 


Equations  II-69,  11-70,  and  11-71  are  not  in  general  satisfied  by  a  solu¬ 
tion  of  the  type  u  =  /  t  +  j  which  satisfies  the  elastic  wave  equation 
(equation  II-7).  However,  solutions  for  equations  II-69,  11-70,  and  11-71 
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of  the  elastic  wave 


can  be  obtained  from  the  solution  u  =  / ^  ^ 

equation  using  the  "correspondence  principle"  (14)  provided  the  boundary 
conditions  for  the  two  problems  are  the  same.  The  "correspondence 
principle"  states  that  the  solution  to  a  visco-elastic  problem  may  be 
obtained  from  the  solution  for  the  corresponding  problem  for  an  elastic 
case  by  applying  the  one-sided  Fourier  transform  to  the  elastic  solution, 
replacing  the  elastic  constants  by  the  corresponding  visco-elastic  moduli 
or  compliances  and  finally  inverting  the  transform.  The  complex  visco¬ 
elastic  compliances  for  the  models  in  fig.  II-7  are  given  as  follows 
(see  Bland,  p.  114-115) 


J(i<n)  =  M  o  p  “  o"  iajTp-p  (Kelvin-Voigt)  11-72 

m  +  anr  m  +  a nr 


J(ia>)  =  I  -  hi  (Maxwell)  11-73 


J(ia)) 


i  “l  i  (,  *4  A 

Mo  +  +  ^  ^  + 


(Standard  Linear  Solid)  11-74 


where  i  =  V^r  and  oo  is  the  frequency  (rad/sec)  of  the  oscillating 
force . 

The  one- 
defined  by 

oo 

f(a>)  =  f  e~iCDt  f(t  ±§)  dt  11-75 

0 


sided  Fourier  transform  of  a  function  such  as  f 


e^) 


IS 


The  inverse  transform  is 
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4  ±  I)  =  I E 


11-76 


eiaJfc  /(cd)  to 


where  R  denotes  the  real  part  of  the  expression  in  brackets.  Therefore, 
the  first  step  in  the  solution  of  a  linear  visco-elastic  dynamic  problem 
is  to  find  the  solution  for  the  elastic  case.  As  an  example  consider  the 
solution  11-10  given  for  a  semi-infinite  elastic  medium. 


II-77 


The  form  of  f  is  determined  by  the  conditions  on  the  boundary  (z  =  0). 

If  an  impulse  of  magnitude  I  is  applied  on  the  boundary  at  time  t  =  0  , 
the  boundary  condition  is 


0(0, t)  =  I  5(t) 


11-78 


where  6(t)  is  the  Dirac  delta  function.  From  equation  11-77  (d) 


|  f’(t)  =  I  5(t) 


II-79 


Integrating  with  respect  to  time 
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II-80 


»-#(*■  &  ® 1  H(t 
H) 
M) 

M) 


-§) 


u  =  §IS 


M 


0  =  I  6( 

where  H  (•-§)  is  the  Heaviside  step  function. 

Let  us  find  the  stress  for  a  visco-elastic  medium  under  the  same 
boundary  load.  Applying  the  "correspondence  principle,"  the  one-sided 
Fourier  transform  of  stress  a  is  found  first 


oo 


a(a>)  =  1  J  e-ia3t  b(t  - 
0 


dt 


II-81 


or 


a(oo)  =  Ie 


-ia^z/C 


11-82 


Next  the  elastic  compliance  is  replaced  by  the  complex  compliance 
j(iao)  of  the  visco-elastic  medium.  Since  C  =*J M/p  equation  11-82 
becomes 


-icoz  pj(ico) 


c(cd)  =  Ie 


The  stress  itself  is  found  from  the  inverse  of  equation  11-83 


11-83 


o(z,  t)  = 


J~  gicD  ^t  -  ^/pj(ico) 

_0 


dec 


11-84 
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For  a  Maxwell  material,  j(ico)  =  j|  -  —  (equation  11-73).  After  substitu¬ 
ting  in  equation  11-84  and  integrating 


_  Mt  r 

o(z,t)  =  Ie  "  §)  + 


where  1^  is  the  first  order  Bessel  function  of  imaginary  argument.  The 
impulse  I  is  propagated  with  velocity  C  and  attenuation  aJ pm/2tj  . 
Similar  expressions  can  be  obtained  for  Kelvin-Voigt  and  standard  linear 
material  by  substituting  their  corresponding  compliances,  J(ico)  ,  in 
equation  11-84  and  performing  the  integrations.  Thus,  the  solution  to 
any  visco-elastic  problem  can  be  found  from  the  solution  of  the  corre¬ 
sponding  elastic  problem  by  using  the  "correspondence  principle"  as  was 
demonstrated  for  the  case  of  the  impulsive  load.  In  Chapter  IV  this 
principle  is  used  for  obtaining  the  solution  for  a  standard  linear  material 
from  the  corresponding  solution  for  the  linear  hysteretic  material. 

F.  Shock  wave  propagation 

In  Section  D  it  was  shown  that  for  a  nonlinear  material  each  stress 


level  Op  travels  at  a  velocity  which  is  given  by 


11-86 


If  the  stress-strain  relation  for  the  material  is  such  that  increasing 

stresses  are  accompanied  by  increasing  values  of  (fig.  II-8),  then 

P 

higher  stress  waves  will  travel  at  higher  velocities  and  eventually  catch 


up  with  the  lower  stress  waves  which  preceded  them  (fig.  II-8).  This  is 


stress  at  z  =  0 


ie 


C 

s 


1/2 


r 

Fig.  II-8 . 


Typical  stress-strain  curve  causing  ’’shocking  up” 
of  the  stress  wave. 
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the  process  of  "shocking  up"  which  is  governed  by  the  fundamental  shock- 
wave  equations  known  as  the  Rankine-Hugoniot  equations  derived  from  the 
equations  of  conservation  of  mass  and  momentum. 

Consider  a  semi-infinite  column  of  material,  having  a  stress-strain 
curve  shown  in  fig.  II -8,  subjected  to  an  instantaneous  step  pressure  applied 
to  its  surface  (fig.  U-9).  Under  the  action  of  the  applied  pressure  P  , 
a  shock  wave  is  initiated  and  propagates  down  the  column  with  some  velocity 
of  propagation  C  .  At  a  time  t  =  t  ,  the  front  will  reach  a  position 
z  =  C  t  ,  as  shown  in  fig.  II-9.  During  an  infinitesimal  time  dt  ,  the 

S  3# 

shock  front  will  have  moved  to  position  z  =  C  (t  +  dt).  The  infinitesimal 
volume,  ACg  dt  ,  where  A  is  the  cross-sectional  area  of  the  element,  will 
then  obtain  a  particle  velocity  u  .  From  Newton's  second  law,  the  relation 
for  impulse -momentum  across  the  shock  front  is  expressed  as 

P  A  dt  =  C  dt  Ap.u 
s  1 


p  = 


C  p.u 
si 


II-87 


where  p^  is  the  initial  density  of  the  material.  The  strain  e  in  the 

infinitesimal  volume  AC  dt  is  defined  as  the  change  in  displacement  per 

s 

unit  length 


u  dt  u_ 

C  dt  "  C 
s  s 


11-88 


By  combining  equation  11-88  with  equation  II-87,  the  shock  wave  velocity 
can  be  expressed  as 


From  equation  11-89,  it  can  be  concluded  that  the  shock,  velocity  is  a 

p 

function  of  the  secant  modulus  —  of  the  material  as  shown  in  fig.  II-8. 
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Stress  Wave  Propagation  Through 
A  Laterally  Constrained  Column  of 
Linear  Hysteretic  Material 

A.  Description  of  model 

The  linear  hysteretic  model  was  first  used  by  Salvadori,  Skalak ,  and 
Weidlinger  (l6)  to  approximate  the  behavior  of  soils  in  a  uniaxial  state  of 
strain.  The  model  is  independent  of  the  rate  of  load  application;  hence, 
energy  dissipation  is  due  only  to  the  compaction  (hysteretic)  character¬ 
istics  of  the  model. 

Consider  the  stress-strain  curve  of  fig.  III-l.  The  stress-strain 

curve  is  a  straight  line  OA  on  initial  loading  to  A  .  The  slope  of  this 

line  defines  an  initial  modulus  M  defined  as  the  elastic  constrained 

o 

modulus.  M  is  related  to  Lame’s  constants  X  and  [j,  ,  or  to  Young’s 
o 

modulus  E  and  Poisson's  ratio  v  ,  as  follows: 

-  »•  +  -  (1 1  i)(l  rW  II1-1 

Upon  unloading  the  stress-strain  curve  is  another  straight  line  AB 
which  defines  a  second  modulus  (it  is  assumed  herein  that  Poisson's 

ratio  is  the  same  for  loading  and  unloading).  If  the  material  is  reloaded, 
it  follows  line  BA  to  A  and  then  continues  along  the  initial  loading 
line  AC  .  In  the  limit  when  approaches  Mq  ,  the  medium  will  behave 

as  a  linear  elastic  material. 

If  the  initial  peak  stress,  peak  strain  and  residual  strain  are  re¬ 
spectively  denoted  as  a  ,  e  and  e  (fig.  III-l),  the  initial 
loading  modulus  Mq  is  then  given  by: 
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III-2 


M 


o 


CT 

max 

6 

max 


The  unloading  modulus  is  given  by: 


M 


1 


a 

max 


e  -  e 
max  r 


III-3 


The  propagation  velocity  of  an  initial  loading  stress  wave  is  given  by: 


c0  =  p  III-1* 

where  p  is  the  mass  density  of  the  material.  The  propagation  velocity 
of  an  unloading  or  subsequent  reloading  stress  wave  is  given  by: 

q  =  in-5 

B.  Boundary  load 

The  dynamic  boundary  load  used  with  this  model  for  blast-type  problems 


is  a  pulse  characterized  by  a  peak  stress  at  a  shock  front  and  an  exponen¬ 
tially  decaying  behavior  thereafter  (fig.  III-2)  given  by 


P(t)  =  PQ  e_t/T 


where : 


III -6 


P  =  the  peak  applied  pressure 
o 

t  =  time,  the  independent  variable 

t  =  the  exponential  time  constant  (time  at  which  pressure  has 
decayed  to  O.368  Pq^ 

C.  Formulation  of  the  problem 

Consider  a  plane  compression  wave  of  general  nature  (arbitrary  rise 
and  decay)  propagating  into  a  column  of  linear  hysteretic  material  ex¬ 
tending  from  z  =  0  to  z  =  00  .  The  distribution  of  pressure  along  the 
column  length,  at  a  given  time  t  ,  might  be  as  shown  in  fig.  HI-3-  The 


48 


Overpressure 


P(t) 


z 


Fig.  III-3.  One-dimensional  wave  conditions. 
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pressure  varies  from  zero  at  the  front  of  the  wave,  point  A,  to  a  maximum 
located  some  distance  behind  the  wave  front,  point  B,  to  the  decayed  over¬ 
pressure  existing  at  the  surface,  point  C.  In  the  portion  of  the  wave  AB 
the  pressure  is  continuously  rising  and,  hence,  this  portion  of  the  wave 
will  propagate  with  velocity  Co  corresponding  to  the  loading  branch,  OA, 
of  the  stress-strain  curve  in  fig.  III-l.  The  one -dimensional  wave  equa¬ 
tions  governing  the  motion  in  zone  AB  are  therefore: 


a? =  c2  at2 

o 


in-7 


s2u  _  i  a2u 
Sz2  e2  St2 

o 


III-8 


where  a  is  the  stress  in  the  direction  of  wave  propagation,  u  repre¬ 
sents  the  particle  displacement  in  the  direction  of  wave  propagation, 
and  u  is  the  particle  velocity.  In  this  zone  the  behavior  is  the  same 
as  that  of  a  linear  elastic  material;  the  stress,  particle  velocity, 
particle  displacement,  and  strain  are  related  by: 


a  =  C  pu 
o 

u  =  C  :e 
o 

u  =  CQ/edt  III -9 

where  e  is  the  strain  in  the  direction  of  wave  propagation.  In  the 
portion  of  the  wave  BC,  the  stress  is  less  than  at  B  and  is  continuously 
decreasing.  This  portion  of  the  wave  will  propagate  with  velocity  C 
corresponding  to  the  unloading  branch,  AB,  of  the  stress-strain  curve  in 
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fig.  III-l.  The  governing  equations  in  this  zone  are: 


Sz2  =  C2  ^2 


d2u  1 

Sz2  C2  St2 


III-IO 


III-ll 


The  response  of  the  medium  in  this  zone  is  not  governed  by  elastic 
theory  because  the  unloading  stress-strain  curve  is  offset  from  the 
origin.-  Since  velocity  is  greater  than  Cq  ,  the  unloading  wave 
front  BC  eventually  overtakes  the  loading  wave  front  AB.  The  result¬ 
ing  interaction,  known  as  internal  reflection,  will  cause  a  decrease  in 
both  peak  particle  velocity  and  peak  stress.  This  process  continues 
with  time  as  will  be  shown  in  the  following  section. 

D.  Response  of  the  medium  to  a 

discontinuous  surface  overpressure 

A  semi-infinite  body  of  linear  hysteretic  material,  such  as  that 

in  fig.  III-3,  is  loaded  uniformly  at  its  surface  by  a  pressure  pulse 

described  by  P(t)  .  Surface  displacement  and  stress  are  assumed  to  be 

zero  before  the  application  of  the  load.  Strains  are  considered  small  and, 

therefore,  density  is  considered  constant.  The  surface  pressure  pulse  is 

suddenly  applied  at  the  time  t  =  0  and  it  remains  constant  at  a  value 

Pq  until  a  time  t  =  t^  ;  the  pressure  then  suddenly  drops  to  a  value 

P^  and  remains  constant  until  time  t  =  co  .  The  more  general  surface 

“*"fc  /t 

pressure  (fig.  III-2),  p(t)  =  PQe  '  ,  will  be  analyzed  in  the  next 
section.  At  the  time  t  =  0  and  at  z  =  0  a  wave  begins  to  propagate 
into  the  body  as  shown  in  fig.  III-4.  For  the  time  0  <  t  <  t^  ,  a  step 
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loading  wave  of  pressure  Pq  propagates  into  the  body  with  a  velocity 

C  .  The  stress  in  zone  0  of  fig.  Ill -4  is  determined  by  the  boundary 
o 

conditions  at  the  surface  and  is  given  as 


a  =  P 
o  o 


III -12 


The  particle  velocity  is  determined  from  the  elastic  theory  (equation 
III -9)  and  is  given  by: 


P 

•  o 

Uo  ~  pC 


III -13 


At  the  time  t  =  t^  ,  a  step  unloading  wave  of  pressure  Pq  -  P^ 
begins  to  propagate  into  the  body  with  a  velocity  which  is  the  slope 

of  the  line  separating  zones  0  and  1  in  fig.  Ill -4.  The  stress  in  zone  1 
is  determined  from  the  boundary  conditions  at  the  surface  and  is  given 
as : 


ai=  P1 


III-14 


The  particle  velocity  in  zone  1,  u^  ,  is  determined  from  an  equation  ex¬ 
pressing  the  uniqueness  of  u^ 


u,  =  u  -  £u  _ 
1  o  ol 


III-15 


where  ZSu^  is  the  change  in  particle  velocity  caused  by  changing  the 


stress  from  a  to  cr_ 
o  1 


^ol  ”  pC 


a  -  ct  P  -  P 
o _ 1  o _  1 


PCn 


III-16 


'1  r  1 

Substituting  equations  III-13  and.  III-16  in  equation  III -15 ,  one  obtains 
for  u^ 
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Ill -17 


*i  -  w  k +  V-*  (p° '  pi} 

o  [_  1 

The  unloading  step  wave  Pq  -  P^  traveling  at  velocity  overtakes  the 
loading  step  wave  Pq  at  the  point  zg  at  the  time  tg  .  From 
fig.  HI -4  it  can  be  seen  that: 


where 


z2  ='=0*2  -  <^(*2  "  V 


Cot2  -  Vs  '  Vl 


t2  -  9>\ 


Cl/Co 

^  =  Cl/Co  -  1 


III-18 


HI-19 


At  time  t  ,  all  the  material  in  zone  1,  from  z  =  0  to  z  =  C  t  , 

d,  O  c. 

is  at  a  stress  of  cr^  and  has  a  particle  velocity  of  u^  .  The  material 
below  z^  is  unstressed  and  at  rest.  It  may  be  considered  that  the 
material  in  zone  1  is  a  body  of  elasticity  striking  the  lower  undis¬ 

turbed  mass  of  elasticity  Mq  .  The  resulting  elastic  interaction  is 
governed  by  the  equations  of  continuity  of  stress  and  particle  velocity 
across  the  boundary.  The  transmitted  wave  travels  downward  with  velocity 
Cq  ,  and  the  reflected  wave,  which  divides  zones  1  and  2,  travels  upward 
with  velocity  .  The  particle  velocity  in  zone  2,  for  the  transmitted 
wave  is 


III-20 


where  is  the  stress  in  zone  2.  The  particle  velocity  for  the  re¬ 

flected  wave  is 
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III-21 


»  •  X\ 

*2  "  "l  “  pC1 


where  c  is  the  reflected  stress.  From  continuity  of  stress  across  the 
K 


boundary: 


a.  +  o'  =  c 
1  R  2 


Substituting  equation  III-22  in  equation  III-21  obtains 


III -22 


u  =  u 
2  1 


III -22 


From  equations  111-17,  III -20,  and  III -23  it  is  found  that: 


’2  ■  P1  +  “(Po  -  V 


111-24 


where : 


CL/C  -  1  , 

l7  o  1 


C../C  +  1 

l7  o 


III -25 


Since  a  >  0  for  all  cases,  the  upward  traveling  wave  is  a  reloading  wave. 
Substituting  equation  111-24  in  equation  III-20  yields  the  particle  veloc¬ 
ity  in  zone  2 


u  =  -~[p_  +  a(  P  -  P, ) 

2  PC  L  1  0  1  J 


111-26 


The  reflected  wave  reaches  the  surface  (z  =  0)  at  time  t^  which  is 


determined  from  fig.  III-4  as 


=  t2  +  (t  -  tp/2 


-  *2  -  h 


III -27 


C.  +  C  ^ 

1  o  t  1 

-  C  i  a 
1  o 
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The  stress 


Conditions  at  time  t^  are  similar  to  those  at  time  t^  . 
is  defined  from  the  boundary  condition  at  t 


°3  =  P1 


111-28 


The  particle  velocity  may  be  written  down  by  analogy  to  u^  ,  replacing 
uq  by  u2  and  o^  by  a 2  ,  respectively,  in  equations  III-15  and 
III-16  obtains 


u  =  u  -  Ziai  „ 
3  2  23 


.  g2 
U2 


HI-29 


Substituting  equations  111-24,  111-26,  and  111-28  in  equation  III-29,  one 


obtains  for  u. 


U  =  -i-  P  +  -  (P  -  P  ) 

3  pC  |_  1  1  +  a  o  l'J 


HI-30 


The  reflection  from  the  moving  wave  front  at  t^  ,  fig.  III-4,  may  be 
analyzed  in  the  same  way  as  the  reflection  at  t^  •  The  unloading  wave 
traveling  at  velocity  C-^  overtakes  the  loading  wave  at  the  point  z^ 
at  the  time  t^  .  From  fig.  Ill -4  it  can  be  seen  that 


z4  ~  c0t4  “  Ci^4  ”  ^ 


*4  -  e 


c  t,  A 

_  +•  -  . 1 _ i.  =  ^  4 • 

Co  3  Cx  -  Co  a  a  *1 


HI-31 


The  particle  velocity  for  the  transmitted  wave  in  zone  4  is 

.  ct4 

u,  =  —z— 

4  PC 


III-32 
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•where  0 ^  is  the  stress  in  zone  4.  The  particle  velocity  for  the  re¬ 
flected  wave  is 


III -33 


From  equations  III-28,  III-30,  III-32,  and  III-33  it  is  found  that 


a4  =  P1  +  a2(P0  -  Px) 


111-34 


Substituting  equation  111-34  in  equation  III -32  for  particle  velocity  in 
zone  4, 


h  *  a2(po  -  V 


HI-35 


The  above  process  can  be  continued  indefinitely  by  finding  the  stress 
at  the  surface  from  the  boundary  condition,  computing  the  new  velocity, 
taking  into  account  the  change  of  stress,  and  then  analyzing  the  inter¬ 
action  from  the  unloading  wave  overtaking  the  initial  loading  wave.  From 
equations  III-18,  III-31,  III-24,  111-34,  III -26,  and  III -35,  one  can 
deduce  the  following  general  expressions  for  arrival  time  t  ,  stress 

0  and  particle  velocity  u  in  any  even  numbered  zone,  that  is, 
n  n 

n  =  2,  4,  6,.... 


t 

n 


111-36 


0 


n 


n/2 
a  ' 


III-37 


u 


n 


111-38 


Similar  expressions  can  be  deduced  for  any  odd  numbered  zone,  that  is, 
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n  =  1,  3)  5,...,  from  equations  III-27,  HI-14,  III-28,  III-17,  and.  III-30 


.  (l-n)/2  . 

t  =  or  "  t. 

n  1 

III -39 

n 

III-40 

XL  =  P,  +  -  (P  -  P,  ) 

n  pC  1  1  +  a  o  1 

r  o  L  -J 

III-41 

It  is  apparent  from  equations  III-40  and  III-37  that  the  stress 

oscillates  from  P^  for  n  odd  to  the  value  given  by  equation  III-37  for 

n  even,  and  approaches  P^  for  large  values  of  t  since  a  <  1  because 

C,  /C  >1  and 
1'  o 

Limit  a11/2  =  0  111-42 

n  — »oo 

Therefore,  for  large  values  of  t  ,  equation  III-37  yields  the  same  re¬ 
sults  as  equation  III-40. 

a  (z  ,  oo  )  =  P  111-43 

The  limiting  particle  velocity  from  both  equations  111-38  and  III-41  is 

u  (z  ,  CO  )  =  III -44 

It  is  now  desirable  to  write  expressions  for  stress  and  particle 
velocity,  as  functions  of  time,  which  are  valid  in  the  zone  0  <  z  <  CQt  . 
Examination  of  equations  III -37  and  III-40,  for  stress  in  even  and  odd 
numbered  zones,  reveals  the  fact  that  the  stress  in  each  new  zone  is  made 
up  of  values  from  the  previous  zone  plus  a  contribution  due  to  the  wave 
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traveling  upwards  or  downwards  at  velocity  .  The  upwards  traveling 
wave  is  a  reloading  wave,  it  has  a  positive  contribution  on  stress.  The 
downwards  traveling  wave  is  an  unloading  wave  and  it  has  a  negative  con¬ 
tribution  on  stress.  For  instance,  the  stress  in  zone  3>  from  equation 
III-40,  is  =  Pj  •  The  stress  in  zone  2,  from  equation  111-37}  is 

a ^  -  Pj_  +  a(P0  -  P^).  Thus,  the  contributing  wave  is  an  unloading  wave 
of  magnitude  a(PQ  -  P^)  separating  zones  2  and  3*  The  stress  in  zone  4, 

from  equation  111-37}  is  cr^  =  P^  +  a  (P  -  P^)  indicating  a  contribution 

-•  ..  _ 

of  a  (PQ  -  P  )  from  the  reloading  wave  which  separates  zone  3  from  4. 

In  general,  stress  in  zone  n  ^  0  is  equal  to  stress  in  zone  n  -  1  plus  a 
contribution  due  to  the  unloading  or  reloading  wave  which  separates  zones 
n  and  n  -  1.  Stress  in  zone  n  -  1  is  also  made  up  of  stress  in  zone  n  -  2 
plus  a  contribution  from  the  unloading  or  reloading  wave  separating  zones 
n  -  2  and  n  -  1.  Therefore,  stress  in  any  zone  n  is  made  up  of  the  sum  of 
all  waves  traveling  downwards  and  upwards  plus  a  contribution  aQ  =  Pq  from 
zone  0.  The  sum  of  all  waves  traveling  upwards,  from  fig.  III-4  and  equa¬ 
tions  III-37  and  III-40  is 


at 


+  a(p  -  p_ )  +  a2(p  -  Pn )  +  or^P  -  p.. ) 

CM  v  o  1'  v  o  1/  K  o  1J 


+  a  (p  - 
o 


pl> 


+  an(p  - 


V 


111-45 


where 


z 


surf ace . 


is  the  travel  time  of  the  reloading  wave  from  point  z  to 
Equation  111-45  can  be  written  in  the  following  form, 

CD 


the 


III-46 
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where 


0  , 


= 


<  V0^ 


n 


an(Po  -P1),t>  t1/an 


111-47 


The  sum  of  all  waves  traveling  downwards  (unloading  waves),  from  fig.  III-4 
and  equations  III-37  and  III-40  is 


Ft  - 


-)  =  -  (Po  -  Pl>  -  a<Po  -  Pl> 

-  a2(Po  -  P1)  -  ...  -  an(Po  -  Px) 


or  using  equation  111-47 


F  ‘  ^ 


oo 


=  -  (P  - 


pi> 


-E 


f  t 
n 


n=l 


(po  -  v  -  g  r  -  c- 


III-48 


111-49 


Therefore,  stress  in  any  zone  n  /  0,  fig.  Ill -4,  can  he  expressed  as 


°(z’v  -ait+-1-Fr-,r  P° 


-  G(t  +  y  -  G(t  -  i^)  -  <Po  -  Pl>  +  Po 


III-50 


P1  +  G|t  +  |-  -Gt-^ 


To  include  zone  0  and  thus  write  an  expression  for  stress  which  is  valid 
in  the  zone  0  <  z  <  CQt  ,  one  must  replace  P^  with  a  function  P(t) 
so  that 


6l 


p(t)  = 


p  ,  t  <  t, 
o  5  1 


Hi-51 


p,  ,  t  >  t 


1 


This  function  is  the  original  input  wave  P(t)  corrected  by  the  travel 
time  f—  ,  that  is,  Pft  -  ^— ]  .  Therefore 

°1  l  Cl/ 


a(z,t)  =  pft  ‘  +  Gr  +  fM '  G(t  ' 


III-52 


As  an  example  we  will  compute  the  stress  at  some  arbitrary  point  z, 
fig.  III-4,  at  a  time  t  =  t,_  . 


°<Z’V  =  P(t5  -  f-  ♦  G  t5  +  §-  -  G  t5  -  ~~ 


III-53 


p/t,-  -  )  is  the  value  of  the  input  wave  at  ft,.  -  ~ ]  ,  from  fig.  III-4 


5  Cl 


5  Cl 


this  value  is  Pn  .  G(t._  +  —  is  the  sum  of  all  the  reloading  waves 

1  l  ^  Clf 

from  t  =  0  to  |t,_  +  —  j  ,  from  equation  III-46 

CD 


111-54 


To  evaluate  equation  111-54,  one  must  use  equation  111-47  in  the  following 
way 

t., 


for 


11  “  1  ’  *5  +  >  V  ’  t3  ’  fi45  +  £4  '  «<po  -  V 


for  n  =  2  ,  t5  +  ^  =  t5  ,  f2|t5  +  §-j  =  a2(Po  -  Px) 

for  n>3  ,  tj.  +  1—  <  =  t„  ,  fr}(tl.  +  rr~ \  =  f  =0 


5  cl  a3  7’  3  5  ‘  ci»  n 
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therefore 


a(P 


V 


a2(P 


pl) 
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The  last  term  in  equation  III -53 


waves  from  t  =  0  to  | t_  -  — 

5  C, 


,  which  is  the  sum  of  all  the  unloading 
can  be  evaluated  in  the  same  manner 


G 


111-56 


from  equation  111-47 

for  n  =  1  ,  tj. 

for  n  >  2  ,  t<- 

therefore 


f„  t 


215 


0 


G/t  -  =-)  =  a(Po  -  Pl) 

-L  i 


III-57 


Substituting  equations  III-55,  III-57,  and  the  value  of  P|t 
in  equation  III-53>  one  obtains  for  o(z,t^) 


z 

Cl 


c(z,t5)  =  P1  +  cr(Po  -  Px) 


111-58 


which  is  the  stress  in  zone  4  obtained  from  equation  III-37. 

The  particle  velocity  time  history  can  be  derived  in  the  same  manner  by 

considering  the  particle  velocity  in  each  new  zone  of  fig.  III-4  to  be  made 

up  of  values  from  the  previous  zone  plus  a  contribution  f  (equation  111-47) 

due  to  the  wave  traveling  upwards  or  downwards.  Since  C,  >  C  ,  from  the 

1  o  5 
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interaction  equations  the  reflected  (reloading)  particle  velocity  will  have 
the  same  sign  as  the  incident  (unloading)  particle  velocity.  Therefore,  the 


contributions  of  all  the  reloading  waves  on  particle  velocity  are 


-G  |t  + 


III-59 


and  the  contributions  of  all  the  unloading  waves  are 


(*  '  °i  )  P0-Pl 


III-60 


The  contributions  of  the  original  input  wave  should  be  written  in  a  form 
that  the  following  conditions  will  be  satisfied 


Particle  velocity  due 
to  original  input  wave 


pC^  when  t  <  t1 
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P  -  P, 
o  1 


when  t  >  t. 


These  conditions  are  satisfied  if  the  contributions  of  the  original  input 
wave  on  particle  velocity  are  expressed  as 


111  +P  f-i _ i- 

PC±  o  1  pCQ  PCX 


111-62 


as  can  be  checked  by  direct  substitution.  The  velocity  time  history,  from 
equations  III-59,  III-60,  and  III-62,  becomes 


u(2>t)=^  p  ,1r  tr  1  +  P°k-^ 


111-63 


E.  Response  of  the  medium  to  a  con¬ 
tinuously  variable  surface  pressure 


Consider  the  decreasing  surface  pressure  as  shown  in  fig.  III-2.  The 
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pressure  p(t)  rises  from  zero  to  Pq  at  t  =  0  .  For  t  >  0  ,  the 
pressure  decays  continuously  to  zero.  This  surface  pressure  may  be  re¬ 
garded  as  a  step  wave  of  amplitude  Pq  plus  a  series  of  infinitesimal, 
negative  steps,  dP  .  Each  pressure  change  dP  may  be  treated  in  the  man¬ 
ner  of  the  single  pressure  change  from  Pq  to  P^  of  the  previous 
section.  When  the  surface  pressure  shown  in  fig.  Ill -2  is  applied  to  the 
boundary  of  a  semi-infinite  body  of  linear  hysteretic  material,  there  will 
result  a  term  like  P(t  -  z/Cj)  (equation  III-51)  where  in  this  case 
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There  will  also  be  a  series  of  contributions  like  G(t  +  z/Cj)  in 
equation  III-52.  In  the  case  of  the  continuously  decreasing  surface  pres¬ 
sure  each  term  of  the  series  in  equations  III-46  and  III-49  will  itself  be 
a  series.  For  a  time  t  ,  the  sum  of  contributions  of  the  first  term 
f^(t),  from  equation  111-47,  will  be 


fx(t) 


111-65 


Integrating  equation  III-65  over  all  times  for  which  contributions  of  the 
form  f-^(t)  are  different  from  zero,  ■  one  obtains 


f-^t)  =  a 


-  -rr  dt  =  aP  -  ap(at) 
dt  o  '  ' 


III -66 


In  the  same  manner  contributions  of  the  form  f  (t)  may  be  written  as 
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where 


P(ant)  =  PQe  T  III-68 

From  equations  III-46,  III-67,  and  III-68  the  contributions  of  all  the 
reloading  and  unloading  waves,  G(t  +  z/c^),  become 


G 


III-69 


Substituting  equations  III-64  and  III-69  in  equations  III-52  and  III-63  one 
obtains  the  following  expressions  for  stress  and  particle  velocity  for  the 
continuously  variable  surface  pressure . 


a(z,t)  =  Pq 


00 


n=l 


HI-70 


u(z,t) 


P 

O 


pc. 


2  a 


n 


n=l 
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F.  Computer  program 

A  computer  code  is  available  for  the  numerical  evaluation  of  equations 
HI-70  and  HI-71.  The  results  of  the  code  computations  are  valid  for  any 


positive  t  and  for  any  a  in  the  range 

0  <  a  <  1.0  III-67 

The  computations  are  made  for  successive  times  at  selected  depths  so  that 

stress-time  and  particle  velocity-time  wave  forms  can  be  constructed. 

Code  output  is  provided  in  two  forms : 

Form  1  Attenuation  of  peak  vertical  stress  and  particle  velocity 
with  depth 

DEPTH  (ft)  SIGMAX  (psi)  VMAX  (fps) 

Form  2  Stress  and  velocity  time  histories  at  specified  depth 

TIME  (sec)  SIGZZ  (psi)  VELOCITY  (fps) 

where 

SIGMAX  =  Maximum  vertical  stress  at  a  given  depth 
VMAX  =  Maximum  particle  velocity  at  a  given  depth 
SIGZZ..=  Vertical  stress 
VELOCITY  =  Particle  velocity 

The  input  variables  for  the  code  consist  of  the  following: 


Variable  _ Description _  Format. 

EL  Loading  modulus  (psf)  -  M  (fig.  III-l)  E  10.2 

EU  tfriloading  modulus  (psf)  -  M^(fig.  III-l)  E  10.2 

POP  Peak  applied  overpressure  (psi)  -  PQ(fig.  Ill -2)  E  10.2 

DM  Mass  density  (slugs/cu  ft)  -  p  E  10.2 

TO  Exponential  time  constant  (sec)  -  T(fig.  III-2)  E  10.2 

IZ  Maximum  number  of  depth  increments  +  1  E  10.2 

(Continued) 


67 


Variable 

Description 

Format 

DZ 

Length  of  each  depth  increment  (ft) 

E  10.2 

DT 

Time  increment 

E  10.2 

TEND 

Time  after  arrival  of  the  front  at  each  depth 
at  which  the  program  ends  (sec) 

E  10.2 

The  first  five  input  variables  need  no  explanation.  The  selection  of 
and  JZ  specifies  the  uniform  spacing  and  total  number  of  evenly- 
spaced  depths  at  which  output  is  required.  TEND  provides  a  fixed  amount 
of  total  time-history  to  be  computed  at  each  output  depth;  since  it  is 
referenced  to  the  time  of  arrival  of  the  front  at  each  depth,  computed 
durations  at  each  depth  will  be  the  same.  DT  regulates  the  number  of  . 
history  points  to  be  computed  between  the  time  of  arrival  and  TEND  at 
each  output  depth. 

The  program  is  written  in  FORTRAN  II  card  language.  Computer  time  on 
a  GE  225  for  a  typical  problem  is  approximately  three  minutes.  A  com¬ 
plete  listing  of  the  program  and  an  example  of  typical  resulting  output 
is  presented  in  Appendix -I II. 

Figure  III-5  shows  computed  stress-time  histories  at  three  different 
depths  for  an  example  problem.  This  figure  shows  the  principal  features 
of  the  linear  hysteretic  model,  that  is,  attenuation  of  peak  stress  and 
the  lengthening  of  the  duration  of  the  stress  pulse  with  depth.  The 
attenuation  of  peak  particle  velocity  is  equal  to  that  of  stress. 
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Stress  (psi) 


Fig.  III-5.  Typical  stress-time  histories  for  several  depths;  linear 
hysteretic  model. 
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APPENDIX  III 


LINEAR  HYSTERETIC  MODEL 
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o  o  o  o  o  o  o  o  u  L)  u  o  o  o  o  o  o  ojo  O  O  O  CJ  OICJ  OOOOOOOOOCJ 


*F iz*  OCT  9  /OCT 


CARD  *A AU+  FLOATING  POINT 


program  title 

STRESS  WAVE  PROPAGATION  THROUGH 
.LINEAR  HYSTERETIC  MATERIAL 


LATERALLY  CONSTRAINED  COLUMN  OF 


PROGRAM  NO. 
PROGRAMMER 

short  titlf 


41-Z5-093 

rohani,r 

ID  LINEAR  H 


HYSTERETIC  MODEL 


ASSUMED 

STRESS 


STRAiN 
RELAT ION 


loading 


.  UNLOADING 


_ _ STRAIN _ 

OUTPUT  "  ~ 

PEAK  STRESS  AND  PARTICLE  VELOCITY  aT  DIFFERENT  DEPTHS f ATTENUAT I  ON 
OF  PEAK  STRESS  AND  VELOCITY] 

STRESS-TIME  AND  VELOCITY-TIME  HISTORY  AT  DIFFERENT  DEPTH 
REOUIRED  INPUT  DATA  “ 

EL^LUaD I NG  MODULUS ( PSF } _ 

EUsItNLOADlNP  MODULUS  I  PSF  J  “ 

_ POpsPbAK  OVER  PRESSURFIPSI i 

~  DMawASS  DENSITY tSLUGS/CU. FT. ] 

TCsEXPOMENTlAL  time  CONSTANT [ SEC] 

I/sTUTaL  no.  OF  DEPTH  IMCREMENTS+I 
D2  =  DtfPTH  I  NCREMENT  [  FT .  ] 

DTsTIME  INCPEMENTISEC1 

TEND-TIME  after  arrival  of  THE  FRONT [ SEC  1 _ 

glossary 

alpha^the  parameter  which  causes  the  peak  stress  to  attenuate 
CUslJNLOADING  WAVE  VELOC T TY [ FPS ) 

CL=LU*DING  WAVE  VELOC T TY [ FPS ] 

SlGMAXsPEAK  stress  f  PS T ]  ~ 

2=VAKI able  depthift. ] 

S!GZZ  =  STRESS  tPSI )  ~ 

VELOCITY IFPSJ 

VMAX=MAX  VELOCITY [FPS]  ~ - 

TPsNUN DIMENSIONAL  ARRIVAL  T-IME* 


AL  ui-  l  tit  r  HUNT  l  SfcC  J 

tfinal=tend+ta 
end  of  glossary 

READ  10* El, EU, POP, DM, TO, TZ,DZ,DT»TENn 


FORM AT [7bl 0.2] - 

CL  =  SC'RTFIEL/DM] 

CU  =  ^QRTF  l  EU/DM ]  - 


I  c 


c 

c 

c 

c 


ALPHAS  H.-tCL/CUJ  l/ll.+lCL/CU]  ]  45" 

PRINT  12 j  EL,EU,PM,  CL.CU,  ALPHA  _  45 j 

12  format iihi,2ix,5hel  =  fio.i,4h  psF.3y.5HEu  =  fig.i,4h  psf.3x.5hDm  wwt 
Is  F3 » 1» 16H  SLGS.PER  CU.FT»,//»21X,5HcL  =  F10.3.4H  FPS.3X.5HCU  s  FI  45? 
20.3/4H  FPS.3X.8HALPHA  =  F6 . 2, // , 21X, 68HATTENUAT I  ON  OF  PEAK  VfcRTICA  45* 
3L  STRESS  AND  PARTICLE  VELOCITY  WITH  DEPTH* //, 28X , 1QHDEPTH [ FT .), IPX  45* 
4.11HSIGMAX IPSI ] .lOX.lOHVMAy [FPSJ  , //)  45f 

2er'#  46_ 

47 

_  48 


DO  600  1=1,  IZ 
TP=Z/tCL*TQ] 


SIGMAX=POP 


49 

- — - - - - -  50 

Js" 501 

ion  Ks j  +  1  502 

”  TERMSEs-pOp*  tl.- ALPHA  1*1  ALpHA»»J|»ti.-EXPF  I  ( -2 .  *TP»ALPHA'**K  j  7  ti.+A  5i 


1LPHA) ]  I 

IFf  ABSFltERMSEl-U.  i!uU5*PnPl  120,110,11! 


52 


110  continue 


54 


SiGNAXsSiGMAX+TERMSE 
Jr  j  +  1 


GO  TO  100 

120  SIGMAXsSIGMAX+TERMSE 

VMAX=  tSiOMAX*144]  /IDM  +  CL1 


PRINT  ll.Z.SIGMAX.VMAx 


55 

56 
37" 

58 
381 

59 


Zs7+DZ 

11  FORM AT  1 1H  ,26X,F5. 1.15X.F7.2.15X.F7.?] 
60';  continue 


60 

601 

61 

62 

63 

6^ 
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_ lilt _ _ _ 66 

DO  6 U 1  N  =  l,  IZ  - 67- 

TAsZ/CL _ . _  68 

TFTNAL  =  TEi\in  +  TA  - - £9~ 

PRINT  13.Z.TA  693 

130FOphat[1H1,21x,3v3HSTRESS  AMD  VELOCITY  T  1'mE  H I  STORY .  .  //,  22X, - 577 

14HZ  =  F7.2.3HFT..5X.5HTA  =  F7.4.4H  SPC , / / *  02 X , 9HT I  ME  I  SEC  1  * 9X , 1 0 HS I  693 

2GZ7l^SI  J  *9X.1'3hvElOC1tyIfH«!J  ,//] - - - “ - 577 

TgT  A _ _ _  70 

19  SlGZZsPUH*EXPFr-IT-Z/CUJ/TOJ - 7T~ 

VV=IPOP*EXPF[-IT-Z/CUJ/Tnn/IDM*CU]  +  IPOP]/ r DM  +  CL ) -  I  Pop] / IDM  +  Cu J  71 3 

”  VVsW*  144  .  ~  - J2~ 

MsT  _ _ _  73 

2'  TERM=P0P»2.*  lALPHA*»M]»gyPFt  - t  AlPHA**M I »tZtO J * [ bXPP 1  [  alPH'A*»H  ]  +L>  \  74 

1TO+CO] J-EXPF I-I ALPHA++MJ+Z/ [TO+CU] U /2.  75 

0XX=  (  1-2.*  POP]  /  fDM*CU)  3*1  f  ALPHA**M]-.5*[ALPmA**mW[FX!=>F[-IaLPHA#*M)  751 
1 *  T / T  0 ] ]*[EXPF[[ALPHA**M]*Z/IT0*CU1 ] +FXPF r - r ALPHA**M ) *1/ [ TO*CU ] 1 ] ]  75? 

XX  =  XX*144  .  '  '  - - 753^ 

IF!  AbSFlTERN]-O.OOf'5*pnPl  220.210,210  76 

210  CONTINUE  - JJ~ 

SIGZZSSIGZ7  +  TEPM _  7« 

W  =  VV  +  XX  Tiaj- 

_ M  s  M  + 1 _  JQ 

GO  TO  20  qq 

220  SIGZZ=SiGZ7+TEPM  fl1 

VV  =  VV  +  XX  ~  — - gyy 

PRINT  14,T,SIG7Z,VV  82 

14  FORMATtlH  ,21X,F6.4,12X. F7 , 2 , 1 2X , F7 . ?  ]  53“ 

TsT*DT _ _____  833 
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IF [T-TFINAL ]  i 9,19,21  64 

21  Zs?  +  DZ _ ^ _ 85 

601  CONTINUE  ”  86 

_ 30  TO  9 _ 87 

ENO  88 

■  1437,  17766 _ 


»KP  0  0  0  n r  0  0 _ 

/ L'  U  0  0 9  0  00  It  0  r.7 

/ron-tn  n- non 65 _ 

EL  00  0  0'' 0 4 

Eli  0000014 _ 

POP  0 0  u  0  ~ 1 6 

DM  nnon:--2c _ 

TO  0000 '122 

IZ  noon  .24 _ 

DZ  POPP' 26 

DT  P  fl  P  n  '•  3 l _ 

TEND  o ' 00032 

CL  -POOP  34 _ 

SORTF  01.41036  FXT  PROG 

clt  o o  n o ■*  4 _ 

ALPHA  0":  01*042 

1  .  0  0  0  0105 _ 

/LG 012  0  'Til  173 

Z  r>  n  n  p  1 1 4 _ 

0»  000011* 

1  o  n  n  p  1  2 ; _ 

*■'  upi  1  p ' 0 1  1 2 1 

JLf _ ■■oon‘|22 _ 

SIGmax  n  '  n ■ 124 

j  ■"f>rip~t26 _ 

+  000--P  0  '  0  1 2 7 

/“run  p- no  356 _ 

K  I  0  U  0 1.  3  i 

TERMSE  p  fi  1362 _ 

,1*1  ..  100364  EXT  PROG 
EXPF  OT40366  FXT  PROG 

2  .  ^  00  0  37 1: 

AfeSF  P t  4 C  4  4  FXT  F P 0 6 
0 .0005  P  ~  n t 4 :  6 

/  012P  n  ^  n  P  5  ;•  3 _ 

/L011  0  0  "  O’.  474 

V  M  A  X  P  '  P  •.  4 1  u _ 

+  f- o 1 4 4  pno t. 5 ! ; 7 

/■.  DPll  o  *  0 1 5  7 1 _ 

/LP6'P  0  "  0 o  6  r  4 

N  M  n  n  n  5 1 ! : _ 

TA  o  n  o  n  ^l? 

T  F  I  N  A  L  n  -■  n  ;  :  5 1 4 _ 

/ 0  0 1  3  0:00635 

T  nnn5i 6 _ 

/ 1.  0  0 1 9  c  *  n  i.  722 

SIGZ7  n-p.520 _ _ 

VV  oO;'P522 

144.  P-*nl(i22 _ 

M  on  m*24 

/,  Cin2P  P  010  57 _ 

TERM  n* 01026 
XX  n poi-44 
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.5  "001046 

/ 1 1  0  2  2  0  0^01560 
/G0210  0  r  01346 
/ 0 0014  n 001400 
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EL  s  3II C  n  0  0 . 0  PSF  EU  =  l2440000.n  PSF  Dm  =  3.1  FIGS. PER  CU.FT. 

CL  =  1  0  0  0  .  0  Of,  FPS  CU  =  2000.000  FPS  ALPHA  =  0733 

ATTENUATION  OF  PEAK  VERTICAL  stress  and  particle  velocity  with  nERTH 


DEPTM ( FT . ]  SlGMAXtPSlJ  VMAXIFPS) 


0.0 

10  0.00 

4.63 

50.0 

69.91 

3.24 

100 . 0 

5  0*68 

2.35 

150.0 

38.18 

1  .77 

2  0  0.0 

29.87 

1.38 

250.0 

24.20 

1.12 

300.0 

20.22 

0.94 

350.0 

17.32 

0.80 

400.0 

15.15 

0.70 
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STRESS  AMD 

VELOCITY  TIME  HISTORY 

• 

Z  =  0.0 

"FT  .  TT~Z  0.0  000 

T  I  ME [ SEC  1 

SIGZZIPSII  ! 

VELOCITY [FPS] 

0 . 0  U  0  o 

100.00 

4.63 

0.0100 

81.87 

4.11 

0  .  U  2  0  o 

67.03 

3.67 

0  .  0300 

54.88 

3.31 

0.0400 

44.93 

2.99 

0 .0500 

36.79 

2.73 

0 . 0600 

30.12 

2.50 

0.0700 

24.66 

2.31 

0 . 0800 

20.19 

2.14 

0.0900 

16.53 

2.00 

0.1000 

13.53 

1.88 

0.1100 

11.08  ' 

1.77 

0.1200 

9.07 

1.66 

0.1300 

7.43 

1.59 

0.1400 

6.08 

1.52 

0.1500 

4.98 

1.45 

0.1600 

4.0  8 

1.40 

0 . 1700 

3.34 

1.35 

0 . 1 6  0  0 

2.73 

1.30 

0.1900 

2.24 

1.26 

0 . 2  0  U  0 

1.83 

1.22 

0.2100 

1.50 

1.19 

0 . 22  0  0 

1.23 

1.16 

0.2300 

1 . 0 1 

17T3 - 

0.2400 

0.6? 

1.10 

0.2500 

0.67 

1.08 

9.2600 

0.55 

1.06 

0.2700 

0.45 

1.04 

0.2800 

0.37 

1.02 

0 . 2900 

C  .30 

1.00 

0.3000 

0.25 

0.99 

0 . 3100 

0.20 

0.97 

0.3200 

0.17 

0.96 

0.3300 

0.14 

0795 - 

0.3400 

o.ll 

0.93 

0.3500 

0.09 

0792 

0.3600 

0  .  (‘7 

0.91 

0  *  3  /  0  0 

0.16 

0.90 

0.380C 

0  .  f‘5 

0.90 

0 . 3900 

0.04 

0.89 

0.4000 

0.03 

0.88 
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STRESS  AMD 

VELOCITY  TIME  HISTORY. 

Z  =  5G.00FT,  TA  =  0.0500  SEC 

T I  ME  £  SEC  1 

SlGZZtPSlJ 

velocity IFPS] 

-v 

0.05C0 

69.90 

3.24 

0.0600 

58.37 

2.90 

0.0700 

^8.86 

2.61 

0.0800 

41.01 

- 2737 - 

0 . 090c 

34.53 

2.16 

0.1000 

29.18 

1.08 

0.110C 

24.74 

1.82 

0.1200 

21.06 

1.69 

0.13G0 

18.00 

1.57 

0.1400 

15.46 

1.46 

0.1500 

13.34 

1.37 

0.1600 

11.56 

1.29 

0.1700 

1 0  .  0  8 

1.22 

0.1800 

6 .83 

1.15 

0.1900 

7.78 

1.10 

0.2000 

6 .89 

1 .04 

0.2100 

6.13 

1.00 

0.2200 

5.49 

0.95 

0.2300 

4.94 

0.91 

0.2400 

4.46 

0.98 

0.2500 

4 . 0  6 

0.84 

0.2600 

3.7  0 

0.81 

0.2700 

3.30 

0.78 

C . 2800 

3.12 

0.75 

0.29  0  0 

2.89 

0.73 

0.3000 

2.68 

0 .70 

0 . 310C 

2.49 

0.68 

0.3200 

2.32 

0.66 

0.3300 

2.18 

0.64 

0.3400 

2 .  G  4 

0.62 

0.3500 

1.92 

0.61 

0.3600 

1 . 61 

0.59 

0.370Q 

1.71. 

0.57 

0.380  0 

1.62 

0.56 

0 . 39  G  0 

1.53 

0.54 

0 . 4  0  0  0 

1.45 

0.53 

0.410  0 

1.38 

0.52 

0.4200 

1 .32 

0.50 

0 . 4300 

1 . 25 

0.49 

0 . 4400 

1.20 

0.48 

0.4500 

1.14 

0.47 
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STRESS  AMD  VELOCITY  TIME  HISTORY. 

Z  *  1 <c . 0 D F T .  Ti  =  0.1000  SEC 

T I  HE  I  SEC  1  SIGZZIPSI J  7  VELOC I T Y [ FPS J 


0.10  0  0 _ 50. 6R  2.35 


0.1100 

43.21 

- 27T3 - 

0.12G0 

37.01 

1.94 

0.1300 

31.84 

1776 

0.1400 

27.52 

1.64 

0.1500 

23.91 

1.52 

0.1600 

20.88 

1.42 

0.1700 

16.33 

- T772 - 

0.1800 

16.17 

1.24 

0.1900 

14.35 

1.17 

0.2000 

12.80 

1.11 

0.2100 

11.47 

IT oS - 

0.2200 

10.34 

1.00 

0.2300 

9.36 

0 .95 

0.2400 

8.52 

0.91 

0.2500 

7.78 

0.87 

0.2600 

7.14 

0.84 

0.2700 

6.58 

0.80 

0.2800 

6.09 

0.77 

0.290C 

5.65 

0.75 

0.300C 

5.26 

0.72 

0.310C 

4.91 

0.70 

0.3200 

4.59 

0.67 

0.3300 

4 . 3l 

0.65 

0.3400 

4.05 

0.63 

0.3500 

3.82 

0.62 

0.3600 

3.61 

0.60 

0.3700 

3.41 

0.58 

0.3800 

3.23 

0.57 

0.3900 

0T55 

0 . 4  U  0  0 

0.54 

O.41G0 

TTn 

- CHT2 - - - 

0.4200 

2.64 

0.51 

0.4300 

2 .52 

0.50 

C  .4400 

2.4n 

0.49 

0.450  0 

2.29 

OTTS 

0.4600 

2.19 

0.46 

0.4700 

2.10 

- PT45 - 

0.4800 

2.01 

0.44 

0.4900 

1 .93 

0.44 

0.5000 

1.85 

0.43 
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STRESS  AMD 

VELOCITY  TIME  HISTORY. 

Z  =  15C.O0FT.  TA  =  0.1500  SEC 

T I  ME  t  SEC  I 

SIGZZIPSI J 

VELOClTY[pPSJ 

r. 

0.1500 

38.17 

1.77 

0.1600 

33.24 

1.62 

0.1700 

?9 . 1 0 

1.50 

0.1800 

25.61 

1.39 

0.1900 

22.66 

1.30 

0.2000 

20.15 

1.21 

0.2100 

18.02 

1.14 

0.2200 

16.19 

1.08 

0.2300 

14.62 

1.02 

0.2400 

13.27 

0.97 

0.2500 

12.10 

0.92 

0.2600 

11. 04 

0.88 

0.2700 

10.19 

0.84 

G.28UC 

9.40 

0.81 

0.2900 

8.71 

0.78 

0.3000 

8.0  9 

0.75 

0.3100 

7.55 

0.72 

0.3200 

7.05 

0.70 

0.3300 

6.61 

0.67 

0.3400 

0.65  “ 

0.3500 

5.85 

0.63 

0.3600 

5.5p 

0.61 

0.3700 

5.21 

0.59 

0.3800 

4.93 

0.58 

0.3900 

4.68 

0.56 

0.4000 

4.44 

0.55 

0.4100 

4.22 

0.53 

0.4200 

4.02 

0.52 

0.43UC 

3.83 

0.51 

0.4400 

3.65 

0.49 

0.4500 

3.49 

0.48 

0.46CG 

3.34 

0.4  7 

0.4700 

3 .19 

0.46 

0 .4800 

3.06 

C  .  45 

0 . 490C 

2.93 

0.44 

0.5000 

2.81 

0.43 

0.5100 

2.70 

0.42 

0 .5200 

2.59 

0.41 

0.5300 

2.49 

0.41 

0.5400 

2.40 

0.40 

0.5500 

2.31 

0.39 
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STRESS  AND  VELOCITY  TINE  HISTORY. 


Z  =  200  .  OnFT  . 

TT" =  ij.2U00  SEC  - 

T I  ME  t  SEC  1 

SIGZZ 1  PS  T 1  1  VELOCITY [TPS] 

0 .2000 

29.87 

1.38 

0.2100 

26.53  “ 

T72V 

0.2200 

23.69 

1.20 

0.2300 

21.26 

l7T3 - 

0.2400 

19. 18 

1.06 

0 . 250C 

17.39 

1.00 

0.2600 

15.83 

0.95 

0.2700 

14.49 

Q  >90  — 

0.2800 

13.31 

0.86 

0.2900 

12.27 

0782 

0.3000 

11.36 

0.79 

0.3100 

10.55 

0.76 

0.3200 

9.83 

0.73 

0.330C 

9.19 

0.70 

0.3400 

8.61 

0 .66 

D . 350  C 

8.08 

0.66 

0 .3600 

7.61 

0.63 

U  .  37  ti 0 

7.18 

0.61 

0.3800 

6.78 

0.60 

0.3900 

6.4p 

- 73E - 

0.4000 

6.09 

0.56 

0.4100 

5.78 

075 5  ' 

0.4200 

5.50 

0.53 

0 .4300 

5.23 

0.52 

U  .  4400 

4.99 

0.51 

0.4500 

4.76 

g-74-9 

0.4600 

4.54 

0.48 

0.4700 

4.34 

0.47 

0.4800 

4.16 

0.46 

0.4900 

3.98 

0.45 

0.5000 

3 . 82 

0.44 

0 .5100 

3 . 66 

0.43 

0.5200 

3.5? 

0.42 

0.5300 

3  *  3  ft 

0.41 

0.5400 

3.25 

0.40 

0.550  0 

3.13 

0.40 

0 . 5  b  0  0 

-3.01 

0.39 

0.5700 

2.90 

“  07  31 

0.5800 

2 . 6  0 

0.37 

0 .5900 

2.70 

0.37 

0 . 6.0  0  0 

2. 6f] 

0.36 
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STRESS  AMD  VELOCITY  TIME  HISTORY 


Z  =  250 . 0  OFT . 

TA  =  0.2500  SEC 

T I  ME  t  SEC  1 

SIGZZIPSIJ 

VELOCITYIFPSJ 

0  *  2  5  C  (j 

24.20 

1.12 

0.2600 

21.87 

1.05 

0.2700 

19.86 

0.99 

0.2800 

18.12 

0.94 

0.2900 

16.60 

0.89 

0.3000 

15.28 

0 . 85 

0.3100 

14.11 

0.81 

0.3200 

13.08 

0.78 

0.3300 

12.17 

0.74 

0.3400 

11.36 

0.72 

0 . 3500 

10.62 

0.69 

0.36  0  0 

9.97 

0.66 

0.37  0  0 

9.37 

0.64 

0.3800 

8.83 

0.62 

0.3900 

8.34 

0.60 

0 . 4  u  0  0 

7.89 

0.58 

0.4100 

7.48 

0.57 

0.4200 

7.10 

0.55 

0.4300 

6.74 

0.54 

0.440C 

6.42 

0.52 

0.4500 

6.12 

0.51 

0.4600 

5 . 84 

0.49 

0.47  0  0 

5.57 

0.48 

0.4800 

5 . 33 

0.47 

0.4900 

5.10 

0.46 

0.5000 

4.88 

0 . 45 

0.5100 

4.68 

0.44 

0.52  0  0 

4.49 

0.43 

0.5300 

4.31 

0.42 

0.5400 

4.14 

0.41 

0.5500 

3.98 

0.40 

0.5600 

3.83 

0.40 

0.57  0  0 

3.69 

0.39 

0.5800 

3.56 

0.38 

0.5900 

3.43 

0.37 

0 . 6  0  U  0 

3.31 

0.37 

0.6100 

3.19 

0.36 

0.6200 

3.09 

0.35 

0.6300 

2.98 

0.35 

0.6400 

2.88 

0.34 

0.6500 

2.79 

0.34 
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STRESS  AND  VELOCITY  TIME  HISTORY 


Z  =  300.0nFT. 

TA  =  0.3000  SEC 

T I  ME [ SEC  1 

SIGZZ ( PSD  T  VELOCITY! fPSI 

0.3000 

20.21 

0.94 

C . 31U0 

38.53 

0.3200 

17.07 

0.84 

0 . 3300 

15.78 

0.80 

0.3400 

14.64 

0.77 

0.3500 

13.63 

0.74 

0 . 3600 

12.7? 

0.71 

™~~  o.37uo 

11. $1 

- EHTB - 

0.3800 

11.18 

0.66. 

0.3900 

10.52 

0.63 

0.4000 

0.61 

0.4100 

9.38 

0.59 

0.4200 

8.68 

0.57 

0.4300 

8.42 

- lTT^S - 

0.4400 

8.00 

0.54 

0 . 4500 

7.61 

0.53  ' 

0 . 46U0 

7.24 

0.51 

0 . 4  7  0  n 

6.9l 

0.50 

0.4800 

6.60 

0.49 

0.4900 

6.30 

0.47 

0.5000 

6.03 

0.46 

0.510C 

5.77 

0.45 

0.5200 

5.53 

G  .  44 

0.5300 

5.31 

0.43 

0 . 5400 

5.10 

0.42 

0.5500 

4 .90 

0.41 

0.5600 

4 . 7l 

0.41 

0.5700 

4.53 

0.40  “ 

0 . 5800 

4.36 

0.39 

0.5900 

4 . 2  n 

0.38 

0.6000 

4.05 

0.37 

0.6100 

3.91 

UT37 - 

0.6200 

3.78 

0.36 

0.6300 

3 . 65 

0 . 35  ~ 

0.6400 

3.53 

0.35 

0.6500 

3.41 

07T4 - 

0.6600 

3  •  3  0 

0.34 

0.6700  ~ 

3T19  " 

0.33 

0 . 6800 

3.0  9 

0.33 

0.6900 

3.00  . 

0.32 

0.7000 

2.91 

0.32 

TTPICAL  CODE  OUTPUT,  LINEAR  HYSTERETIC  MODEL 


STRESS  AMD  VELOCITY  TIME  HISTORY 


Z  =  350 . 0  OFT . 

.  TA  =  0.3500  SEC 

T I  ME ( SEC ] 

SIGZ7IPST  I  VELOCITYfFPSI 

0.3500 

17.32 

0.80 

0 . 3600 

16.07 

0.77 

0.3700 

14.96 

0.73 

0.38GC 

13.98 

0.70 

0.39CC 

13.00 

0.68 

0.4000 

12.29 

0.65 

0.410C 

11.5? 

0.63 

0.4200 

10.92 

0.61 

0.4300 

10.32 

0.59 

0.4400 

9.78 

0.57 

0.4500 

9.28 

0.55 

0.4600 

8.62 

0.54 

0.4700 

8.39 

0.52 

0.4600 

8 . 00 

0.51 

0 .4900 

7.63 

0.49 

0 . 5  u  0  0 

7.29 

0.48 

0.5100 

6.97 

0.47 

0.5200 

6.67 

0.46 

0.5300 

6.39 

0.45 

0.5400 

6.13 

0.44 

0.5500 

5.88 

0.43 

0.5600 

5.65 

0.42 

0.5700 

5.43 

0.41 

0.5800 

5.23 

0.40 

0.5900 

5.03 

0.39 

0 . 6  u  0  0 

4.65 

0.38 

0.6100 

4.67 

0.38 

0.6200 

4 . 5l 

0.37 

0.6300 

4.35 

0.36 

0 . 6400 

4.20 

0.36 

0.6500 

4.06 

0.35 

0.6600 

3.93 

0.34 

0.6700 

3 . 8  n 

0.34 

- 

0.680  0 

3 .68 

0.33 

0.6900 

3.57 

0.33 

0.7000 

3.46 

0.32 

0.7100 

3.35 

0.32 

0 .7200 

3.25. 

0.31 

0 . 73C  0 

3.16 

0.31 

0 . 7400 

3.0  7 

0.7500 

2.98 

TYPICAL  CODE  OUTPUT,  LIHEAR  HYSTEKETIC  MODEL 


84 


STRESS  AND  VELOCITY  TIME  HISTORY 


Z  =  4  !"■  0 . 0  0  F  T  . 

TA  =  0.4000  SEC 

T  I  ME l SEC  1 

SIRZZIPSIJ  ! 

VELOCITY [FPS] 

0.4u00 

15.15 

0.70 

0.4100 

14.19 

0 .6? 

0.4200 

13.33 

0.65 

0.4300 

12.55 

0.62 

0.4400 

11.64 

0.60 

0.4500 

11.20 

0.58 

0.4600 

10.61 

0.56 

0.4700 

10.07 

”  0.55 

0.4800 

9.57 

0.53 

0.4900 

9.11 

•  0.52 

0.5000 

8.69 

0.50 

0.5100 

8.29 

0.49 

0.5200 

7.9? 

0.48 

0.5300 

7.58 

0.46 

0.5400 

7.26 

0.45 

0.5500 

6.96 

0.44 

0.5600 

6.68 

0.43 

0.5700 

6.41 

0.42 

0.5800 

6.16 

0.41 

0.5900 

5.9? 

- OTTO - 

0 . 6  U  0  0 

5.70 

0.40 

0.6100 

5.49 

0.39 

0.6200 

5.30 

0.38 

0 . 6  3  G  0 

5.11 

0.37 

0.6400 

4.93 

0.37 

0.6500 

4.76 

0.36 

0.6600 

4.60 

0.35 

0.6700 

4.45 

0.35 

0.6600 

4.30 

0.34 

0.6900 

4.17 

0.33 

0.7000 

4 .04 

0.33 

0.7100 

3.91 

0.32 

0.7200 

3.79 

0.32 

0.7300 

3  •  6  8 

0.31 

0.7400 

3.57 

0.31 

0.7500 

3.47 

0.30 

0.7600 

3.37 

0.30 

0.7700 

3.27 

0.29 

0.7600 

3.18 

0.29 

0.7900 

3.10 

0.29 

0.8000 

3.02 
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Stress  Wave  Propagation  Through,  a  Laterally  Confined 
Column  of  Visco-Elastic  Hysteretic  Material 

A.  Introduction 

Dynamic  compression  and  wave  propagation  tests  have  indicated  that 
some  soils,  like  clay,  exhibit  both  strain-rate  and  hysteretic  effects. 
Increase  in  modulus  with  strain-rate  and  lag  of  maximum  strain  behind 
maximum  stress,  an  indication  of  viscous  behavior,  have  been  observed  by 
many  investigators  (17,24)  during  dynamic  compression  tests  of  some  soils. 
Smoothing  of  the  stress  pulse  during  one-dimensional  wave  propagation 
tests  also  supports  the  viscous  behavior  of  such  materials.  In  order  to 
predict  wave  propagation  phenomena  in  such  soils  a  mathematical  model 
should  be  constructed  to  account  for  both  the  hysteretic  energy  loss 
and  the  viscous  behavior  of  soils.  The  viscous  behavior  of  materials  is 
frequently  represented  in  terms  of  rheological  models  consisting  of  linear 
springs  and  dashpots.  Three  types  of  such  models  were  discussed  in 
Chapter  II  and  the  stress-strain-time  relationships  for  each  model  were 
derived  in  terms  of  the  spring  constants  and  the  viscosities  of  the  dashpots 
(equations  11-60,  11-64,  11-68).  It  is  a  difficult  task  to  decide 
which  of  these  models  would  best  describe  the  viscous  behavior  of  soils. 

In  fact,  no  single  model  can  describe  viscous  phenomena  in  all 
soils.  The  limited  research  performed  in  this  area  (19)  has  indicated 
that  the  three  element  standard  linear  visco-elastic  model  (fig.  II-7)  may 
be  used  to  describe  the  viscous  behavior  of  some  clays  and  sands  within  cer¬ 
tain  ranges  of  boundary  and  loading  conditions.  Further  experimental  work 
is  required  in  this  area  if  rational  rheological  models  are  to  be  constructed 
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which  will  describe  the  viscous  phenomena  of  soils  more  accurately  and  over 
a  broader  range  of  boundary  and  loading  conditions.  For  the  present  anal¬ 
ysis  however,  the  standard  linear  visco-elastic  model  will  be  used  to  de¬ 
scribe  strain-rate  effects.  The  "correspondence  principle,"  described 
in  Chapter  II,  will  be  utilized  to  obtain  visco-elastic  solutions  from 
equations  III-70  and  III-71  for  stress  and  particle  velocity  in  a  linear 
hysteretic  half space  subjected  to  a  uniform  time -dependent  pressure  wave 
at  its  surface.  The  solution  for  any  other  linear  rheological  model  may 
be  obtained  by  the  same  procedure.  The  linear  hysteretic  model  is  used  here 
to  account  for  the  hysteretic  energy  loss  of  the  soil. 

B.  Formulation  and  solution  of  the  problem 

The  linear  hysteretic  model  was  analyzed  in  Chapter  III.  The  fol¬ 
lowing  expressions  were  obtained  for  particle  velocity  and  stress  within 
the  medium  (subjected  to  an  exponentially  decaying  input  shock)  as 
functions  of  time: 


u(z,t)  =  Pq 
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where 

u(z,t)  =  particle  velocity  at  depth  z  and  time  t 
o(z,t)  =  stress  at  depth  z  and  time  t 


where 


e  =  2.7183... 

Pq  =  peak  applied  stress 

=  unloading  reloading  wave  velocity 

CQ  =  loading  wave  velocity 
CL  /C  -  1 

a  =  J4  0 

c1/co  +  1 

t  =  exponential  decay  constant  of  the  input  stress  pulse 
P  =  mass  density  of  the  medium 

The  front  of  the  wave  propagates  with  the  loading  wave  velocity  CQ  and 
reaches  the  point  z  at  the  time  t  =  z/Cq  .  The  medium  at  z  is  undis¬ 
turbed  for  times  t  <  z/C  therefore,  the  one-sided  Fourier  transforms 

'  o 

(equation  11-75)  of  equations  IV-1  and  IV-2  are  non-zero  when 


z/Cq  ^  t  <  00 


IV-3 


The  first  step  in  the  "correspondence  principle"  is  to  evaluate  the  one¬ 
sided  Fourier  transforms  of  equations  17-1  and  17-2  within  the  limits 
given  by  equation  17-3 


u(cd) 


c(cd) 
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-icot 
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u(zjt)  dt 


J  e  ia>t  cr(z,t)  dt 

z/C 
'  o 


17-4 
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where  u(^)  and  a(u>)  are  the  one-sided  Fourier  transforms  of  particle 
velocity  and  stress  respectively.  Substituting  for  u(z,t)  and  a(z,t) 
from  equations  17-1  and  17-2  and  integrating  one  obtains 
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Fig.  IV-1.  Rate-dependent  model  for  the  linear-hysteretic  medium. 
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The  second  step  is  replace  the  elastic  constants  Cq  and  in  equations 

17-6  and  IV-7  by  |pJo(io>)J _1/2  and  jjoJ^ico)  ~1/2  respectively,  where 
JQ(ico)  and  J^(ico)  are  the  visco-elastic  complex  compliances  of  the  segments 
OA  and  AB  of  the  hysteretic  stress-strain  curve  shown  in  fig.  17-1.  An 
assumption  will  be  made  here  to  reduce  the  number  of  parameters  and  simplify 
the  analysis,  that  is 
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The  complex  compliance  of  the  standard  linear  model  was  given  in  Chapter  II 
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If  one  introduces  the  notation 
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Equation  IV-l4  can  he  written  in  terms  of  its  real  and  imaginary  parts  as 
follows 

^pJo(io>)  =1^(0))  -  ik2(a>)  IV-15 

where  k^(<u)  and  k2(o)  are  the  real  and  imaginary  parts  of  ypj^(in)  . 

From  equations  IV-11  and  IV-15 ,  the  elastic  constants  C  and  Cn  should 

o  1 


then  be  replaced  by 


IV-16 
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respectively.  Therefore  the  parameters  XQ  ,  X^  ,  X^  ,  Sq  and  in 

equation  IV-8  become 

-i“  ^(aiJ-izkg^J 
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So  =  T  P  [V®5  "  ik2(CDi] 

P 

si  =  F-  (1  ‘  p)  [ki(Q)  '  ik2(“}] 


The  third  and  the  final  step  in  obtaining  the  visco-elastic  solution  is  to 
take  the  inverse  transforms  of  equations  IV-6  and  IV-7  (see  equation  11-76) 
with  parameters  Sq  ,  ,  Xq  ,  X^  ,  and  X 2  given  by  equation  IV -18  and 

then  choose  the  real  parts  of  the  transforms,  that  is 

u(z,t)  =  |  R 

I  oo 


CO 

/ 


elCDt  u(o>)  d P 


17-19 


a(z,t) 


o<a»)  d<n 


17-20 


where  R  denotes  the  real  part  of  the  expression  in  brackets.  From  the 


theory  of  integration  of  a  complex -valued  function  (15)  equations  IY-19 
and  IV-20  can  be  -written  in  the  following  forms. 
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u(z,t)  =  |  J  u(CD)]  dCD 
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00 

R[elCDt  c(©)] 


IV-22 


The  functions  n[eia>t  u(a>)]  and  R  [eia>fc  c(^)]  are  determined  from  equa- 

i  0 

tions  IV-6,  IV-7,  and  IV-18  by  using  Euler's  formula,  e  =  cos  9 
+  i  sin  0  ,  and  selecting  the  real  part  of  each  expression.  The  process 
is  lengthy  but  straight  forward;  as  an  example,  the  real  part  of  the  first 


term  in  equation  IV-6  will  be  calculated  here 
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Substituting  for  Sq  and  5^  from  equation  IV-18  and  simplifying,  one 


obtains 
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where 


A  =  cnzk2(co)  +  35^(05)  (1  -  p) 
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TV-26 


Substituting  Euler's  formula 


ioffi 

e  =  cos  a®  +  i  sin  a® 
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into  equation  IV-24  and  simplifying,  obtains 
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or 
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Following  the  same  procedure  and  substituting  the  values  of 
and  R  [eiujfc  cr(to)]  in  equations  IV-21  and  IV-22,  one  obtains 
ing  general  expressions  for  the  particle  velocity  and  stress 
elastic  hystereti'c  medium  shown  in  fig.  IV-1. 
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and  A  and  B  are  given  by  equations  IV-25  and  IV-26  respectively. 

C.  Evaluation  of  k^(u>)  and  k2(<*>) 

To  evaluate  equations  IV-30  and  IV-31}  one  must  know  k.  (cu)  and  kg  (o^ ) , 


the  real  and  imaginary  parts  of  Jq  ( in )  (equation  IV-lU) .  This  can  be 

accomplished  if  equation  IV-14  is  written  in  the  polar  form  of  a  complex 

number.  The  polar  form  of  a  complex  number,  z  =  x  +  iy  ,  is  given  as 
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z  =  r  (cos  9  +  i  sin  9)  where  r  =  yx  +  y  and  9  =  tan  —  .  Therefore, 
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6  =  tan 
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Since  (cos  0  - 


i  sin  0)1/2  =  e“ie/2  =  cos  ^  -  i  sin  ^  ,  equation  IV-37 


becomes 
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The  real  and  imaginary  parts  of  ‘yp Jq ( i  )  are  therefore 


1^(03)  =  ^P/kQ  "Vr  cos  | 
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k2(03)  =  ]jp/kQ  Vr  sin  | 
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respectively. 

D.  Determination  of  visco-elastic  parameters  for  application  to 
soil  problems. 

There  are  five  material  constants  for  the  linear  hysteretic  visco¬ 
elastic  model  under  consideration  which  have  to  be  determined  experi¬ 


mentally,  they  are 
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The  visco -elastic  parameters  $  and  p,  can  be  determined  from  a  dynamic 
uniaxial  strain  test  in  which  a  step  loading  is  applied  and  the  resultant 
strain -time  history  is  measured  (fig.  IV-2).  The  strain  response  of  the 
three-element  visco-elastic  material  to  a  step  pulse  Aog  ,  from  equation 
11-68,  is  given. .as 
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The  solution  of  this  differential  equation  yields  the  following  expression 
for  strain-time  history.  r 


Ae(t)  = 


IV-43 


Substituting  for  Mo/k-L  and  M^/n  from  equation  IV-4l  and  simplifying 
one  obtains 
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Therefore  M  is  the  secant  modulus  at  the  stress  level  Act  at  t  =  Q 
O  s 

(fig.  IV-2).  The  slope  of  equation  IV-44  at  t  =  0  is  given  as 
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For  large  values  of  time,  equation  IV-44  reduces  to 
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IV-47 


tos  o 


Therefore,  the  visco-elastic  constants  M  ,  p,  ,  and  can  he  obtained 
from  the  measured  strain-time  history  using  equations  IV-45,  IV-46,  and 
IV-47  (See  fig.  IV-2).  However,  due  to  the  fact  that  an  instantaneous 
step  load  cannot  be  applied  in  a  uniaxial  strain  device,  approximations 
(20)  to  laboratory  data  (measured  strain-time  history)  are  needed  for  the 
measurement  of  the  visco-elastic  constants. 

The  parameter  g  in  equation  IV-4l  can  be  expressed  in  terms  of  a 
as 


P  = 


1  -  ot 
1  +  a 


IV-48 


The  compacting  dissipative  parameter  a  ,  is  given  by  the  ratio  of  the 
slopes  at  zero  strain-rate,  that  is,  by  the  ratio  of  the  moduli  Mq  and 
^  (fig.  IV-1) . 
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If  there  were  no  viscous  effects,  the  moduli  Mq  and  would  be  pro¬ 

portional  to  the  loading  and  unloading  slopes  of  the  static  stress -strain 
curve.  The  rate -dependent  character  of  the  material  however,  will  modify 
the  proportionality  to  some  extent  and  will  alter  the  definition  of  a 
(see  equation  III-25).  Seaman  (17)  has  developed  an  iterative  procedure 
for  the  approximate  measurement  of  a  from  the  dynamic  uniaxial  test 
data.  Typical  values  of  ot  obtained  by  Seaman  (17)  are  0.12  for  dry 
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sand  (110  pcf  dry  density)  and  0.02  for  kaolinite  (102  pcf  dry  density 
and  33  percent  water  content)  under  100  psi  dynamic  stress. 
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V 


Stress  Wave  Attenuation  Through  a  One-Dimensional 
Column  of  Nonlinear  Locking  Material 


A.  Description  of  locking  media 

The  dynamic  behavior  of  certain  soils  and  rocks  under  uniaxial  strain 
conditions  may  be  approximated  by  the  stress-strain  curve  of  fig.  V-l. 

This  type  of  material  response  is  referred  to  as  locking  behavior.  The 
loading  branch  of  the  stress-strain  curve  is  defined  by  the  relation 


d2g 

de2 


>  0 


During  unloading  the  stress -strain  curve  is  defined  by 


V-l 
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.  .max 


V-2 


If  the  material  is  reloaded,  it  will  follow  the  vertical  branch  of 
equation  V-2  to  the  previous  peak  stress  level,  a  (corresponding  to 

ITlclX 

e  ) ,  and  from  there  on  the  curve  of  equation  V-l. 
max  5 

Experimentally  obtained  loading  stress -strain  curves  are  often  rea¬ 
sonably  fitted  by  a  parabolic  relation  (21)  of  the  form 


-(*)" 


V-3 


in  which  k  and  n  are  constant  characteristics  of  the  medium.  Because 

d2 

of  the  condition  — ^  >  0  ,  it  is  required  that  k  >  0  and  n  >  1  .  The 

de 

following  treats  only  the  relation  given  by  equation  V-3;  other  analytical 
fits  may  be  dealt  with  in  a  similar  manner. 

B.  Boundary  load 

The  dynamic  boundary  load  considered  is  the  pulse  characterized  by  a 


102 


peak  stress  at  a  shock  front  and  an  exponentially  decaying  behavior  there¬ 
after  (fig.  III-2)  given  by 

P(t)  =  Poe“t//,T  V-4 

•where 

Pq  =  the  peak  applied  pressure 
t  =  the  independent  variable  time 

t  =  the  exponential  time  constant,  time  at  which  pressure  has 
decayed  to  O.368  Pq 
C.  Formulation  of  the  problem 

Consider  a  semi- infinite  locking  medium  (fig.  V-2)  that  is  subjected 
to  a  suddenly  applied  pressure  pulse  given  by  equation  V-4  on  its  entire 
free  boundary.  As  a  result  of  the  pulse,  a  shock  is  propagated  into  the 
medium,  compacting  it  to  a  density,  p  (note  that  p  will  vary  with  depth) 
Due  to  the  nature  of  the  stress-strain  relation  (fig.  V-l),  the  density 
change  will  be  permanent.  When  the  shock  front  has  progressed  to  a  depth, 
z  ,  at  a  time,  t  ,  the  mass  mobilized  (involved)  in  the  momentum  transfer 
(assuming  unit  area)  is 


where  pi  is  the  initial  mass  density  of  the  medium.  Because  of  the 
monotonically  decaying  character  of  the  input  pulse,  the  material  behind 
the  shock  front  is  continually  unloaded  along  the  vertical  branch 
e  =  emax  of  the  stress-strain  curve  (no  strain  recovery).  Consequently, 
the  mobilized  mass,  p^z  ,  acts  as  a  rigid  body  causing  all  the  material 
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Fig.  V-2.  Shock-wave  propagation  conditions  in  locking  media 
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behind  the  shock  front  to  move  with  an  identical  velocity  u(t)  ,  while  the 


free  boundary  is  displaced  by  the  amount  u(t)  .  In  accordance  with  the 


Newton's  law,  the  integrated  impulse,  I 


t 

I  =  J  Poe“t//T  dt  V-6 

0 


applied  in  the  same  period  of  time  equates  to  the  sum  of  the  product  of 
elements  of  the  mobilized  mass  and  their  individual  particle  velocities, 
that  is 


t 

p .  zu  =  f  P  e  dt 

1  J  ° 

0 


V-7 


The  stress,  <?  ,  at  the  shock  front  is  given  by  (see  equation  11-86) 

a  =  p.uC  V-8 

ri  s 


where  Cg  is  the  shock-front  velocity  given  by  (see  equation  11-88) 


V-9 


and,  e  ,  is  the  strain  corresponding  to  a  at  the  shock  front  which 
is  related  to  the  particle  velocity,  u  ,  by, (see  equation  II-14) 


V-10 


Substituting  equation  V-10  into 
arranging,  yields 


z 


dz 

dt 


equation  V-7,  using  equation  V-9  and  re- 


V-ll 


106 


or 


dz  TFo  /-  ~ t /t > 

2  at  -  T7  (1  '  e  > 

rl 

Substituting  equation  V-9  and  its  first  integral  into  equation  V-12 
results  in  the  following  nonlinear  differential  equation 

*  1/2 


V-12 


( 


as 


^  /  (I)  « - 
0 


V-13 


Differentiating  both  sides  of  equation  V-13  with  respect  to  time  yields 


1/2 


(ae) 


-1/2 


de  ,  da 
CT  dt  +  e  dt 


/(f) 

0 


1/2 


dt  +  a  =  P  e 
o 


“t/t 


V-l4 


Since 


de  de  da  da/dt 
dt  da  dt  da/de 


V-15 


Equation  V-l4  becomes 


,  H(  sl/2  1  da 


1  +  -s4s- 

da/de 


/(f) 

0 


1/2 


dt  +  a 


P  e 
o 


-t/r 


V-16 


Using  the  results  of  equation  V-13,  equation  V-l6  becomes 


1/2 


P  T  . 
o  da 

a  dt 


1  + 


a/e 

da/d 


(1 


-  e-V'i 


)  +  a  =  Poe 


-t/T 


V-17 


The  variable  coefficient  in  the  brackets  depends  on  the  stress-strain 
curve.  The  ratio 
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V-18 


£ 

e 


is  the  secant  modulus  and 


V-19 


is  the  tangent  modulus.  The  case 


V-20 


corresponds  to  a  linear  locking  material  (fig.  V-3)* 

From  equation  V-3,  the  stress-strain  relation,  one  obtains 


1  + 


q/e 

dq/de 


V-21 


Substituting  equation  V-21  into  equation  V-17,  it  reduces  to  Bernoulli 
equation  (22). 


*  ¥  i  b  0(1  ■  e't/T) 


+  0  =  P  e 
o 


-t/i 


V-22 


With  the  initial  condition 


q(0)  =  P(0)  =  Pr 


V-23 


the  solution  to  equation  V-22  is 


q  = 


tPq[1  -  e  ^/T] 


2n 

n+1 


n-1 

ir  /  a  -  ^-t/T]n+1 « 


/ 

0 


V-24 
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M 


M. 


M. 


—  —■  Secant  Modulus 
e 

cLq 

—  =  Tangent  Modulus 
ae 

M.  For  Linear  Locking  Material 
z 


Vertical  Strain  (e^) 

Fig.  V-3.  Stress-strain  curve  for  linear  locking  material 


t> 
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Equation  V-24  relates  the  stress,  a  ,  at  the  shock  front  with  the  corres¬ 
ponding  arrival  time.  The  shock- front  velocity,  from  equations  V-3  and 

V-9,  is 


n-1 

2n 


172 


(kp±) 


V-25 


The  shock-front  location,  z  ,  from 'equation  V-12,  is  given  by 


z 


i-  (1  - 


V-26 


Utilizing  equations  V-3 ,  V-24,  and  V-25,  the  shock-front  location  time 
relation  becomes 


Equations  V-24  and  V-27  relate  the  peak  stress,  a  ,  at  the  shock  front, 
to  the  shock-front  location,  z  ,  from  which  one  obtains  the  attenuation  of 
peak  stress  with  depth.  The  corresponding  peak  particle  velocity-depth 
relation,  from  equations  V-8  and  V-25,  is  given  by 


u 


V-28 


D.  Numerical  calculation 

A  numerical  procedure,  based  on  the  conservation  laws  of  mass  and 
momentum  and  described  by  Zaccor  (23),  has  been  adopted  for  the  calculation 
of  peak  stress  and  particle  velocity  attenuation  in  the  nonlinear  locking 
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medium.  This  procedure  treats  soil  models  composed  of  all  nonlinear  load¬ 
ing  branch  forms  for  which  stress-strain  behavior  is  described  by  the  re¬ 


lation 


.2 

H>o 

de 


V-29 


and  a  vertical  unloading  branch  given  by  equation  V-2.  Hence,  the  mate¬ 
rial  defined  by  equation  V-3  is  included  as  a  special  case.  The  applied 
boundary  load  can  be  any  pulse  characterized  by  a  peak  stress  at  a  shock 
front  and  a  constant  or  steadily  decaying  behavior  thereafter. 

The  procedure  has  been  validated  (23)  by  comparison  with  the  closed- 
form  solutions  presented  in  section  C  of  this  chapter  for  a  parabolic  load¬ 
ing  stress -strain  curve  and  an  exponentially  decaying  input  pulse. 

E.  Computational  procedure 

The  computational  procedure  described  by  Zaccor  involves  the  impulse - 

momentum  law  given  by  equation  V-7  and  the  concept  of  the  mobilized  mass 

acting  as  a  rigid  body  behind  the  shock  front.  Material  properties 

(loading  branch  of  the  stress-strain  curve  and  the  initial  mass  density) 

and  boundary  conditions  determine  the  velocity,  C  ,  of  the  shock  front 

s 

and,  hence,  its  location,  z  ,  at  any  time.  The  shock -front  location  at 
any  time  identifies  the  mass,  m  ,  mobilized  in  the  momentum  transfer  up 
to  that  time.  The  integrated  impulse,  I  ,  applied  in  the  same  period 
equates  to  the  total  momentum  transferred,  from  which  one  obtains  the  peak 
particle  velocity  and  stress  at  the  shock  front,  that  is 


mu  =  I 
I 

u  =  — 
m 

a  =  p.C  u 
s 


Ill 


The  step-by-step  procedure  used  in  the  program  is  given  below  to 

illustrate  some  of  the  concepts  involved  in  the  theory. 

1.  Compute  the  shock-front  velocity,  C  ,  at  a  discrete  number  of 

s 

stress  levels  (the  program  will  accept  50  values  of  C  ),  from  the  loading 

s 

branch  of  the  stress-strain  curve  (for  0  <  a  <  Pq)  using  the  relation 


C  = 
s 


1/2 


Pi6 
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and  store  the  information  in  tabular  form. 

2.  Take  a  small  time  increment  At  . 

3.  Compute  the  impulse,  l(l)  ,  at  the  time,  t  =  At  ,  as  the  product 
of  the  time  increment  At  ,  and  the  average  value  of  the  overpressure 
pulse  in  that  time  increment 


1(1)  =  At  TP  +  P(At)]/2 


V-32 


where  P(At)  is  the  value  of  overpressure  pulse  at  t  =  At  . 

4.  From  the  shock-front  velocity  versus  stress  table  of  step  1  deter¬ 
mine  the  average  shock-front  velocity 


C 

s 


C  (P  )  +  C  (At) 
s  'o'  s 

2 


in  the  interval  At  as  if  the  stress  in  the  medium  at  the 
were  the  same  as  the  surface  overpressure  at  the  time  t  = 
is  the  shock -front  velocity  corresponding  to  o  =  Pq  and 
shock-front  velocity  corresponding  to  a  =  P(At)  . 

5.  Compute  the  shock-front  location  ,  at  the  end 
increment  At 


V-33 

time  t  =  At 

At  .  C  (P  ) 
s  o' 

C  (A  )  is  the 

S  X 

of  the  time 
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At 


v-34 


6.  From  l(l)  obtained  in  step  3  compute  a  value  for  the  particle 
velocity  u^  at  the  end  of  the  time  increment  At  (see  equation  V-30) 


1(1) 
u,  =  s  ' 
i 


V-35 


where  p^z^  m°bilized  mass  at  the  time  t  =  At  . 

7.  From  u,  obtained  in  step  6  and  C  (At)  from  step  4,  compute 
1  s 

the  stress  at  the  shock  front  at  the  time  t  =  At  and  depth  z^ 

(see  equation  V-30) 


G1  "  PiCs(At)^l 


V-36 


This  value  of  stress  is  now  a  more  accurate  representation  of  the  true 
stress  in  the  medium  at  z^  than  the  assumption  made  in  step  4,  that  is, 
a  =  P(At)  . 


8. 


Compute  a  new  value  of  Cg  using 


C 

s 


Cs(Po)  +  Cs (ol} 


V-37 


where  C  (cm )  is  the  shock-front  velocity  corresponding  to  the  current 

S  JL 

stress  a,  found  in  step  7.  C  (on)  may  be  determined  by  interpolation 
1  s  J- 

in  the  table  computed  for  step  1. 

9.  Start  again  at  step  5  and  continue  to  repeat  this  iterative 
procedure  until  all  values  agree  (within  a  specified  convergence  criterion) 
with  those  on  the  previous  iteration  and  store  the  values  obtained  on  the 
final  iteration,  that  is,  z1  ,  .  Note  that  and  ct  are  the 

peak  particle  velocity  and  stress  at  depth  z^  . 
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10.  Take  the  next  time  increment.  The  total  elapsed  time  is  now 
t  =  2At  . 

11.  Compute  the  impulse  for  the  second  time  increment  as 

1(2)  =  At[P(At)  +  P(2At)]/2  V-38 

where  P(2At)  is  the  value  of  overpressure  pulse  at  t  =  2At  . 

12.  Compute  the  total  impulse  l(2At)  ,  to  the  time  t  =  2At  as 

l(2At)  =  1(1)  +  1(2)  V-39 

13*  Assume  an  average  shock-front  velocity 

C =  C  (c  )  v-4o 

s  si 

where  C  (a  )  is  the  shock-front  velocity  corresponding  to  the  final  value 

S  _L 

of  as  found  in  step  9,  and  compute  the  shock-front  location  at 

the  end  of  the  total  elapsed  time  t  =  2At 


z 


2 


z,  +  C  At 
1  s 


V-4l 


14.  From  l(2At)  obtained  in  step  12  and  z^  obtained  in  step  13, 
compute  a  value  for  the  particle  velocity  u^  at  the  time  t  =  2At  and 
location  zg 


V-42 


15.  Repeat  the  outlined  analogous  procedures  starting  at  step  7 
until  values  of  z^  ,  u^  ,  are  obtained  which  satisfy  the  convergence 
criterion  in  step  9,  and  then  index  to  the  next  time  increment  starting  in 
step  10. 
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F.  Computer  program. 


The  above  procedure  is  programmed  in  FORTRAN  II  language  for  a  GE  225 

computer  with  8K  memory. 

The  input  variables  for  the  computer  code  (program)  consist  of 

the  following: 

Variable  Description  Comments 

ROE 

Initial  mass  density 
(slugs/cu  ft)  - 

Singular  parameter 

DT 

Time  increment  (sec) 

-  At 

Singular  parameter 

sig(i) 

Stress  (psi)  -  a 

Discrete  stress  values 
from  stress-strain  curve 

EPS  (i) 

Strain  (in. /in.)  -  e 

Discrete  strain  values 
from  stress -strain 

curve 

T 

Time  (sec)  -  t 

Discrete  values;  t  =  0 

at  P  =  P 
o 

P 

Overpressure  (psi) 

-  P(t) 

Discrete  values;  P  =  P 
at  t  =  0  ° 

The  input  sequency  and  formats  are  as  follows: 

_ _ Input  Card _  _ Format 


(A) 

(ROE,  DT) 

(F7.2, 

4x, 

F7.5) 

(B) 

[SIG(I),  EPS(I),  1=1,  50] 

(F7.2, 

5X, 

F8.6) 

(c) 

(T,P) 

(F8.6, 

5X, 

F7.2) 

Card  form  (b)  is  repeated  50  times  (for  the  50  discrete  points  on  the 
stress-strain  curve).  Card  form  (C)  is  repeated  (pulse  duration/DT  +  l) 
times,  i.e.,  discrete  values  of  the  input  pressure  pulse  are  read  in  and 
stored  in  tabular  form.  The  code  conveniently  accepts  these  values  at 
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times  which  are  continuous  multiples  of  the  time  increment  DT  .  Note  that 
SIG(l)  =  EPS (l)  =  0.0  (see  fig.  V-l);  SIG(50)  =  PQ  ,  the  peak  overpres¬ 
sure:  and  T  =  0.0  when  P  =  P  . 

o 

The  output  variables  are  given  in  two  groups  as  follows: 

Group  1  Shock-front  velocity  calculation 

STRESS  STRAIN  SHK-VEL 

(psi)  (in. /in.)  (fps) 

Format  (F7.2)  Format  (F8.6)  Format  (F7.2) 

Group  2  Attenuation  of  peak  particle  velocity  and  stress  with  depth 

Z(ft)  V(fps)  STRESS  (psi) 

Format  (f6.2)  Format  (f8.3)  Format  (F9.3) 

0.  Additional  program  information 

The  program  contains  a  function  subprogram,  VELOC,  that  inter¬ 
polates  linearly  in  the  shock-front  velocity  versus  stress  table  calcu¬ 
lated  in  step  1.  The  following  error  indications  are  given  by  sub¬ 
program  VELOC: 

a.  STRESS  TABLE  OVERFLOW 

This  message  is  printed  out  if  the  maximum  value  of  the  overpressure 
pulse  exceeds  the  peak  stress  input,  i.e.,  Pq  >  SIG(50)  . 

b.  STRESS  DATA  OUT  OF  ORDER 

This  message  is  printed  out  if  the  stress -strain  input  are  not  in  an 
ascending  order,  i.e.,  SIG(l)-SIG(l-l)  <  0  . 

H.  Run  time  information 

Computer  time  for  a  typical  problem  is  about  15  minutes.  A  complete 
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FORTRAN  II  listing  of  the  program  and  typical  input  and  the  resulting 
output  for  an  example  problem  is  presented  in  Appendix  V. 
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APPENDIX  V 


LOCKING  MODEL 
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o  olo  olo  o  o  o  o  olo  o  o  o  o  o  o  o  o  o  o  o  o  n  o  o  o  olo 


TITL 

ONE-DIMENSIONAL  STRESS  WAVE  PROPAGATION  IN  LOCKING  AND  DISSIPATIVE 
MED  I A 

PROGRAM  NO.  41-Z5-111  DATE  AUGUST  1967 

Programmer  rohanTTb  Tls 


ASSUMED 

STRESS- 

STRAIN 

relation 


UNIAXIAL 

STRAIN 

DATA 


STRAIN 

GLOSSARY  - 

R0E  =  MASS  DENSITY  OF  SO  I L [ SLUGS/CU . FT  1  _ 

Sl0m=STR6SS  AT  POINT  I  ON  STRESS-STRAIN  CURVEIPSJJ 
EPSII 1«STRAIN  AT  POINT  I  ON  STRESS-STRAIN  CURVEt IN/IN] 
UIIJsSHOCK  FRONT  VELOCITY  AT  STRESS  LEVEL  SIGIIHFPSJ 

_ DT*T I  ME  I NCREMENT  [ SEC ] 

P* OVER PRESSURE  tPSIJ 
T*TIME  A  VAR  I  ABLE  t  SEC ] 


PULSsIMPULSISEC-PSIJ  - 

_ Zf DEPTH  A  VARIABLEIFT. ) 

UBARaAVE.  SHOCK-FRONT  VELOC I TY t FPS 1 
VsPARTICLE  VELOCITYIFPS1  { MAX ] 

STRESS=STRESStMAX)  IPSll” 

VELOC=DUMMY 


INPUT 


ROE,DT,SIGl I] ,EPSII1  ,T,P 
OUTPUT _ 

MAX.  STRESS,  VELOCITY  AS  A  FUNCTION  OF  DEPTH  [“ATTENUATION] 


COMMON  SIG(50],EPS[50)»U(50] 


1  READ  lu 0 , ROE, DT 
LOO  F0RMAT[F7.2,4X,F7.“5J 
PRINT  99, ROE, DT 

990FORMATllHl,l2X,33HSHOCK-FfiONT  VELOCITY  CALCULAT I  ON . » /, 15X, 6HR0E  * 
2F7.2,16H  SLGS  PER  CU.FT. ,1X,5HDT  =  F7.5,6H  SEC.  , /, 15X, 6HSTRESS* 6X 
3,6HSTRAIN,4X»7HSHK  VEL, / , 16X, 5H l PS  I] , 6X, 7H [I N/ IN ] , 3X , 5H [ FPS ] » // J 
DO  1C  I=l,*r 


READ  101, SIGH], EPSII]  - 

Util  eSQRTF I t SIG  t I  1 *144 . } / t (EPS  t I  I ♦ . 0Q0Q011 »R0E1 1 


59 


10  PRINT  l02,SIGm,EPSUJ#um 
101  FQRMAT[F7.2,5X,F8.61 


102  FORMAT l 1H  , 1 4X, F7 , 2, 5X. F8 . 6 , 2X, F7 . 2 J 
PRINT  88 


61 

62 


880FORMAT[1H1,12X,54HATTENUATION  OF  PEAK  PARTICLE  VELOCITY 
_ STRESS. ,/,l6X,5HZtFTl,7X,6HV[FPSl, IPX, HHSTRESStPSIl,! 


63 

641 


AND  PEAK  S  642 
_ 643 


Z  =  0. 

STRESSES IG  t  50 1 


V« tt STRESS J *144. T/ t l ROE UU 1501 ] 
PRINT  134, Z, V, STRESS 


104  F0RMAT[15X,F6.2,6X»F8.3»9X,F9.3] 
READ  103,T,P 


103  F0RMAT(F6.6,5x,F7.21 

PULSsq, 


644 

645 


646 

647 


648 

65 


ZL  =  0 
PL=P 


READ  1  '  3  ,  T  »  P 
PLL  =  P 


66 

67 


68 

70 


71 

111. 


UUsVELOCIPLLl 

^0_PULSsDT* I PL  +  P 1  * . 5+PULS 

UBARs l VELOC  t  PL  1 ♦ VELOC l P 1 i*.5 


SSsC. 


35  Z=UBAR*DT+ZL 

Vs[144.*[PULS/tR0E*Z 1 ] 1 


761 

73 


74 

741 


75 

76 


STRESSsR0E*V*UU*ll./144. 1 
IFtABSFISTRESS-SSl-.OOll  50,50,40 
40  SSsSTRESS 

UBARsVELOCtSSl 


UUsVELOCtSSl' 


77 

78 


781 

79 


790 

791 


GO  TO  35 

50  PR  INT  104,  Z,V ,  STRE SS_ 
UB AR* VELOC (STRESS 
ZL  =  Z 


READ  103, T,p 

PULSs[pT»(PLL*Pl 1/2.+PULS 


80 

81 


811 

821 


83 

84 


PLLsp 
GO  TO  35 


END 

FUNCTION  VELOCtSlSMAl 


PARAMETER  DECK  LOADED  BY  SUBROUTINE 
COMMON  Sl6(5Ul,PPSt501,Uf501 


841 

86 


87 

86 


85 

87 


SsSIGMA 


DO  2  1=2,50 
IFtSIGtI)-SI  2,3,4 


2  CONTINUE 


68 

89 


90 

91 


92 

93 


PRINT  106 
CALL  EXIT 
108  FORMAT  l 


23H  STRESS  TABLE  OVERFLOWTT 
4  Tl=  S I G I  1 1-SIG( I  — 1 1 


IFITU  55,55,70 
55  PRINT  105 


To 5  format 


94 

95 


96 

99 


100 

101 


26H  STRESS"  DATA  OUT  OF  ORDER.  1 


102 

103 


7  •  T2 


um 


Ut 1-1  1 


104 
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TYPICAL  INPUT 
CAED  A 


.44 


I 

ROE 


O.Ol 

I 

DT 


1 1 1 1 1 1 1 1 1 1 1 1 f  f  m  t , , , , , , , , , , , , , , , , , , , , i7iT777TT?T7T7iT777777T77?77TTTT;7?777777? 
22222222222222222222222222222222222222222222222222222222222222222222222222222222 
33H333333333I333333333333333333333333333333333333333333333333333333333333333333 
4444||4444444444444444444444444444444444444444444444444444«444444444444444444444 
55555555555555555555555555555555555555555555555555555555555555555555555555555555 
86B666S66G666GG66G6GBGE666GG6G66GBG6G6GGGBG66BE68G6BB6GG66G6BG6G6B66GGBG6666G6G6 
77777777777777777777777777777777777777777777777777777777777777777777777777777777 

»n|B.il.l.«,|,n,ni„l|„,lllgl|I,|,,1|||>>|8I!8|||18||gi|  8III,  llllllmil  m, 
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TYPICAL  INPUT 
CAED  B 


,  0.00  0.00  N 

I  I 

SIG  (1)  EPS  (1) 

0  0 10 II  0  0  0  0  0  0  1  0  II  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  8  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 

I  !  3  <  !  t  7  «  9  ID  n  C  13  H  13  B  17  II  B  M  7T  72737*  »»  V  7!!3  W  It  7733  34  333S  37  313940  <1  «««4546  47  4I49SC5I  H5JWS5W  57  W  HKfil  67  OW65KS7  St  (9  70  71  72  73  74  75  71  77  71  7SK)  f 

m  t  m  i  n  mi  i  ii  i  m  n  1 1 1  m  i  m  1 1  ii  1 1  m  it  i  m  i  h  1 1 1 1 1 1 1  n  i  m  n  1 1 1 1 1 1  u  1 1 1  m 
22222222222222222222222222222222222222222222222222222222222222222222222222222222 
333|333333333|333333333333333333333333333333333333333333333333333333333333333333 
44444444444444444444444444444444444444444444444444444444444444444444444444444444 
55555555555555555555555555555555555555555555555555555555555555555555555555555555 
666SG66SEEB6R6GE6CGG6E66666G666GE666666E6G6EE6CEE6BEB6S6$6E6G6G66SG66666BEB6SG66 
77777777777777777777777777777777777777777777777777777777777777777777777777777777  | 

0  8  0  1  B  8  8  8  8  8  8  8  8  1  8  8  0  !  8  8  8  B  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  0  8  8  8  8  8  8  8  8  8  0  0  8  8  8  8  8  8  8  8  8  B  8  8  8  8  8  8  8  8  8  !  8  8  !  8  0  0  8  0  It  I  | 
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VI 


Stress  Wave  Propagation  Through  a  Laterally  Constrained 
Column  of  Nonlinear  Hysteretic  Material 

A.  Description  of  model 

The  nonlinear  hysteretic  model  was  recommended  by  Hendron  (25)  on 
the  basis  of  his  study  of  the  behavior  of  sand  under  uniaxial  strain 
conditions.  The  model  is  rate -independent  and  energy  dissipation  is  only 
due  to  the  compaction  (hysteretic)  characteristics  of  the  model. 

A  stress-strain  curve  for  the  nonlinear  hysteretic  model  is  shown 
in  fig.  VI-1.  The  initial  loading  curve  is  given  by  the  relation 

o'  =  As  vi -1 

The  unloading  curve  is  given  by  the  relation 

o  =  B(e  -  er)n  VI-2 

where  A,  B,  and  n  are  constants  characteristic  of  the  medium  and  e 

r 

is  the  residual  strain.  If  the  material  is  reloaded,  it  will  follow  the 

curve  of  equation  VI-2  to  its  previous  maximum  stress  level,  a  ,  cor- 

5  max  ’ 

responding  to  e  ,  and  from  then  on  the  curve  of  equation  VI-1.  The 
nonlinear  hysteretic  model  is  identical  with  the  linear  hysteretic  model 
(Chapter  III)  when  n  =  1.0. 

The  initial  loading  tangent  modulus  Mq  for  the  nonlinear  hysteretic 
model,  from  equation  VI-1,  is  given  as 

do  H™1 

Mo  =  —  =  Ane  (Tangent  modulus)  VI -3 
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Fig.  VI-1.  Typical  stress- strain  curve 
for  nonlinear  hysteretic  model 


The  unloading /reloading  tangent  modulus  ,  from  equation  VI- -2,  is 
given  by 

I L  =  Bn(e  -  VI-4 

For  the  conditions  of  shock  wave  propagation  the  initial  loading  wave 
velocity  is  determined  by  using  the  secant  modulus  of  the  material  (see 
equation  II-89).  The  initial  loading  secant  modulus  is  given  as 

M  =  —  =  Ae11  ^  (Secant  modulus)  VI-5 

o  e 

B.  Boundary  load 

The  dynamic  boundary  load  used  for  this  model  is  a  pulse  characterized 
by  a  peak  stress  at  the  shock  front  and  an  exponentially  decaying  behavior 
thereafter  (see  fig.  II-2)  given  by 

P(t)  =  P  e_t/T  VI-6 

o 

where 

P  =  the  peak  applied  pressure 
o 

t  =  time,  the  independent  variable 

t  =  the  exponential  time  constant  (time  at  which  the  overpressure 
has  decayed  to  O.368  P  ) 

C.  Solution  of  the  problem 

A  numerical  method  developed  by  Heierli  (26)  will  be  used  for  the 
solution  of  the  nonlinear  hysteretic  model.  This  procedure, often  referred 
to  as  the  method  of  impulses,  assumes  that  (a)  the  applied  pressure  pulse 
is  divided  into  a  finite  number  of  steps  containing  a  certain  amount  of 
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impulse,  (b)  the  dynamic  stress -strain  properties  of  the  material  are 
known  for  all  conditions  of  loading,  i.e.,  initial  loading,  unloading, 
and  subsequent  reloading  and  unloading,  (c)  each  increment  of  stress 
change  propagates  at  a  velocity  consistent  with  the  secant  modulus  for 
the  stress  increment,  (d)  the  impulse  is  conserved  on  an  incremental 
basis  and  continuity  is  maintained.  Using  the  conservation  laws  of  mass 
qnd  momentum,  an  expression  for  the  change  in  particle  velocity  and  stress 
can  be  obtained  which,  in  conjunction  with  the  propagation  velocity-stress 
relation  (obtained  from  the  stress -strain  properties  of  the  material),  can 
be  used  to  determine  the  stress  in  refracted  and  reflected  waves  which 
propagate  away  from  the  intersection  between  two  waves.  It  then  becomes  a 
bookkeeping  process  to  compute  the  position  of  all  the  waves  on  a  space- 
time  diagram  (such  as  the  one  shown  in  fig.  III-4)  as  well  as  the  state 
(i.e.  stress  and  particle  velocity)  existing  in  each  zone  of  the  diagram. 
D.  Particle  velocity-stress  relation 

The  particle  velocity-stress  relation  is  determined  from  the  con¬ 
sideration  of  the  conservation  laws  of  mass  and  momentum  for  a  region  of 
material  subjected  to  a  step  change  Ac  in  stress.  This  step  change 

in  stress,  which  moves  from  position  z  to  z  +  As  in  a  time  At  , 

•  •  • 

causes  an  abrupt  change  in  particle  velocity  from  u  to  u  +  Au  as 
shown  in  fig.  VI -2.  The  strain  also  changes  from  e  to  e  +  Ae  .  Assum¬ 
ing  that  the  position  z  is  attached  to  the  moving  particle  (Eulerian 
coordinate  system),  the  conservation  of  mass  for  the  element  of  fig.  VI -2 
yields 

p  Az  =  (p  +  Ap)  [Az  -  Au  At]  VI -7 
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Surface  overpressure 


Fig.  VI' -2.  Propagation  of  loading  waves;  relationships  for 

step  change  in  stress 
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where  Ap  is  the  change  in  density  due  to  the  step  change  Aa  in  stress. 
Equation  VI-7  can  he  simplified  to  the  following  equation 


Ap.  Az  =  Au  At  VI -8 

P  +  Ap 

From  conser'vation  of  momentum 

p  Az  Au.  =  Aa  At  VI-9 

where  paz  is  the  mass  involved  during  the  momentum  transfer  process  and 
A p  At  is  the  applied  impulse  during  the  time  t  =  At  •  The  absolute 
(Lagrangian)  wave  velocity  C  at  which  the  step  ao  travels  is  given  by 

c  -  i +  “  vl-10 

Using  equation  VI-10  in  equation  VI-9  gives 

Aa  =  p  Au.  (C  -  u)  VI -11 

Combining  equations  VI -8,  VI-10,  and  VI -11  to  eliminate  Au  and  solving 
for  A a  gives 

^  <c  -  ^ 

Tie  step  change  Aa  is  assumed  to  be  related  to  the  strain  change 
Ae  by 


Aa  =  M  Ae  VI-13 

s 

where  M  is  the  constrained  secant  modulus  for  the  stress  increment 
s 

At  as  shown  in  fig.  VI -2.  Since  one -dimensional  Eulerian  strain 
is  used 
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VI  -14 
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loading  stress  pulse  will  propagate  into  the  medium  as  a  shock.  The 
velocity  of  propagation  of  the  initial  loading  stress  wave  is  then  de¬ 
termined  from  the  secant  modulus  given  by  equation  VI-5.  In  terms  of 
stress,  equation  VI-5  can  be  written  as  (using  equation  VI-l) 


M 


s 


n-1 

n 

a 


VI -20 


Substitution  of  equation  VT-20  into  equation  VI-18  results  in  the  follow¬ 
ing  expression  for  the  particle  velocity-stress  relation  for  the  initial 


loading 

n+1 

2n 

£  _  - —  (initial  loading)  VI -21 


where  cr  corresponds  to  the  peak  stress  at  the  wave  front  (CT  =  Pq  at  the 

surface).  During  unloading  from  a  stress  level  a  to  a.  ,  the  incre- 

J  "’-L 

mental  change  in  particle  velocity  can  be  obtained  from  equation  VI-19 
where  Mg  is  determined  from  equation  VI -4.  In  terms  of  stress,  equa¬ 
tion  VI-4  becomes 

1  n-1 

M1  =  Mg  =  nBn  c  n  VI -22 


Substituting  from  equation  VI-22  into  equation  VI-19  and  integrating  from 
aj-l  °j‘  5  results  in  the  following  expression  for  the  change  in 

particle  velocity  due  to  a  step  change  (u.  —  o.  _ )  in  stress  during  un- 
loading  or  reloading 


,0-1 


(unloading/reloading ) 


VI -23 
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E.  Wave  velocity- stress  relation 

The  wave-velocity,  Cq  ,  for  the  initial  stress  wave  (a  shock)  is 
determined  from  the  secant  modulus  given  by  equation  VI-5,  that  is. 


c 

o  \  p1 


n-l\  V2 


VI-24 


In  terms  of  stress,  equation  VI-24  becomes 

( .  ¥\1/2 


c  = 


.  n 

A  a 


(initial  loading) 


0  \ci  / 

where,  as  in  equation  VI-21,  a  corresponds  to  the  peak  stress. 

The  wave  velocity  for  the  unloading  or  reloading  stress  wave  (from 
ck  to  aj_-]_)  is  determined  from  equation  VI-2  (the  unloading/reloading 
stress-strain  relation),  that  is 


VI-25 


'3,3-1 


p-i  -  y 

(€o-i  - V  pi 


vi-26 


Note  that  the  initial  mass  density  is  used  in  equation  VI-26;  this  is 

satisfactory  when  the  strains  are  small.  From  equation  VI-2 

1 
,n 


e .  =  €  + 
3  r 


e .  t  =  e 

<3-1  r 


0 . 

_3 

B 


a  .  , 

.3-1 


B 


VI-27 


Substituting  from  equation  VI-27  into  equation  VI-26  and  simplifying 
results  in  the  following  relation 

1/2 


cJ5j-i 


B 


2n  /  0 .  n  -  0  .  \ 

: /  3-1 _ i__  \ 


1  1 

n  n 

\  CTj-l  " 


VI-28 
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The  radical  in  equation  VI-28  may  he  expanded  (l8)  in  terms  of 

and  a  .  ,  -  a  .  to  give 
j-l  3 


°J-1 


+  a . 
3 


I  g.i-l 
1 
n 

\  ^-1 


+  0 


n-l 

2n 


1 


1  -  i- 

1  -  i 

f  2  .  I] 

(° 3-1  ~  i 

1  12 

n 

r  n) 

va-i +  Oo  j  - 

VI-29 


VJ  .  u  . 

If  — ^ — - — -  <  l/k  the  first  term  of  equation  VI-29  is  sufficient  to  give 

an  accuracy  of  one  percent.  Then  the  wave  velocity,  equation  VI-28,  for 
the  above  conditions  becomes 


n 


2 


n-l 

2n 

( unloading/ reloading ) 


VI-30 


F.  Intersection  relations 

Consider  the  space-time  diagram  shown  in  fig.  VI-3.  The  input  pres¬ 
sure  time  is  approximated  by  a  series  of  step  changes  at  the  appropriate 

times  t,  ,  t0  ,  t0  ,  ....  The  first  step  is  given  as  P  at  time  t 
J.  ^  3  o  o 

followed  by  a  decrease  in  pressure  equal  to  £P  =  Pq  -  P  at  time  t^  . 
The  stress  in  the  first  zone  of  the  space-time  diagram  is  equal  to  Pq 
(the  applied  pressure  during  the  first  time  interval) . 


VI- 31 


The  particle  velocity  in  the  same  zone  is  obtained  from  equation  VI-21 
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Overpressure 


Fig.  VI -3.  'Space -time  diagram 


137 


which  expresses  the  relation  between  stress  and  velocity  across  a  shock-front 


n+1 

p2n 

o 


1 


VI-32 


At  time  t^  a  decrease  in  pressure  occurs  at  the  surface  and  causes  an  un¬ 
loading  wave  to  travel  into  the  medium.  The  stress  behind  this  unloading 
front  (in  zone  2)  is 

a2  =  P1  VI-33 


The  particle  velocity,  u2  ,  in  zone  2  is  determined  by  computing  the  change 
in  velocity  £&  due  to  the  step  change  Pq  -  P^  in  stress,  that  is 


u2  =  Up  +  Au12 


VI-34 


where  au-,  p  is  determined  from  equation  VI-23  by  substituting  equa 

tions  VI-31  and  VI-33  for  a.  and  a.  ,  respectively 

3  J 


l+n  1+n 


VI-35 


The  unloading  stress  wave,  Pq  -  P1  ,  governed  by  the  unloading  stress- 
strain  curve  (equation  VI-2),  travels  with  a  velocity  determined  from  equa¬ 
tion  VI-30  which  is  higher  than  the  initial  loading  wave  velocity  given  by 
equation  VI-25 .  The  unloading  wave  therefore  overtakes  the  loading  wave  and 
the  intersection  of  these  waves  produces  a  reloading  (reflected)  wave  which 
propagates  back  toward  the  surface  and  another  (refracted)  wave  which  con¬ 
tinues  on  to  greater  depths.  The  region  between  these  two  wave-fronts 
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has  a  single  stress  and  a  single  particle  velocity. 

The  intersection  equations  can  be  derived  by  considering  the 
uniqueness  of  particle  velocity  at  the  point  of  intersection.  If 
regions  above  and  below  the  interaction  of  the  two  waves  are  denoted  by 
a  and  b  ,  respectively,  and  the  common  region  by  c  as  shown  in 
fig.  VI-4,  the  particle  velocities  above  and  below  the  common  region 
c  axe  given  by 

u  =  u  -  Au  (above)  VI -36 

c  a  ac 

and 


uc  =  ^  +  All^  (below) 


VI-37 


where  An  is  the  change  in  particle  velocity  caused  by  changing  the 

SIC 

stress  from  cr  to  a  and  All  is  similarly  defined.  These  particle 

9.  C  DC 

velocity  increments  are  determined  from  equation  VI -21  for  initial  loading 
or  equation  VI -23  for  reloading  or  unloading.  The  stress  is  determined 

from  the  simultaneous  solution  of  equations  VI-36  and  VI-37. 


The  stress  and  the  particle  velocity  u^ 


of  the  third  zone 


in  the  space-time  diagram  in  fig.  VI -3  will  be  determined  utilizing  the 
intersection  equations  VI-3 6  and  VI-37;  c2  ^2  ’  given  ^y  equations 

VI-33  and  VI-34,  respectively,  correspond  with  c  and  u  in  fig.  VI -4; 

9/9/ 

•  •  • 

and  u^  correspond  with  oc  and  uq  ,  and  c^  and  u^  are  zero  be¬ 
cause  the  intersection  occurs  at  the  leading  edge  of  the  shock-front.  The 

•  • 

particle  velocity  change  Au  is  given  by  equation  VI-23  and  All  is 

9/C  DC 

given  by  equation  VI -21.  From  equation  VI -21 
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(reflected  wave) 


(unloading  wave) 


<vV 


Fig.  VI-4.  Intersection  of  unloading  wave  with  loading 

wave  in  space  time 


Xl+1 

2n 


VI-38 


From  equation  VT-23 


1  +  n 
2n 


VI-39 


Substituting  from  equations  VI-38  and  VI-39  into  equations  VI-37  and 
VI-36,  respectively,  and  eliminating  u^  =  \  gives 


n+1 


•  p^n  B 


2n 


1  +  n 
2n 


VI-40 


Equation  VI-40  can  be  solved  for  c^  .  The  particle  velocity  u^  can 
then  be  determined  from  equation  VI-38 


n+1 

2n 


VI-41 


The  process  just  demonstrated  can  be  continued  indefinitely;  succes¬ 
sively  finding  the  stress  at  the  surface  from  the  boundary  conditions, 
computing  the  new  velocity,  taking  into  account  the  change  of  stress,  and 
then  analyzing  the  interactions  from  unloading  waves  overtaking  the  initial 
loading  wave .  To  obtain  more  accuracy  the  time  between  tQ  ,  t^  ,  tg  ,  ... 
should  be  divided  into  finer  divisions  as  shown  in  fig.  VI- 5 •  The 


l4l 


Fig.  VI-5*  Space-time  diagram  for  a  fine  mesh 


lk2 


intersections  in  this  case  are  not  restricted  to  the  simple  pattern  of 
fig.  VI- 3-  The  intersection  equations  for  an  intersection  which  is  not 
along  the  leading  shock- front  are  also  based  on  equations  VI-36  and  VI-37. 
The  stress  and  the  particle  velocity  u^  are  no  longer  zero  and 

AU-ac  and  are  both  determined  from  equation  VI-23- 

G.  Computer  program 

A  computer  code,  based  on  the  foregoing  analysis  as  related  to  the 
problem  depicted  in  fig.  VI-5?  is  available  to  compute  the  stress  and 
particle  velocity  histories  at  selected  depths  in  a  homogeneous  non-linear 
hysteretic  medium.  The  code  was  originally  written  at  Stanford  Research 
Institute  in  ALGOL  language  (26)  and  has  been  translated  to  FORTRAN  II 
language  for  a  GE-225  computer  with  8k  memory.  The  mechanics  of  the 
code  are  described  in  references  18  and  26. 

Output  parameters  are  given  in  terns  of  the  following  nondimensional 
quantities 


q  _ 

S~Po 

VI-42 

1-3 

II 

VI- 43 

• 

V  =  T" 

VI-44 

U1 

II 

H  1 

oIn 

VI-45 

o 


1 


1 


VI-46 
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where  S  ,  T  ,  V  ,  and  Z  are  the  nondim en sional  stress,  time,  particle 
velocity,  and  depth  respectively;  is  particle  velocity  corresponding 
to  the  peak  overpressure  given  by  equation  VI-32  and  CQ  is  the  wave 
velocity  of  a  shock  wave  with  a  stress  cr  =  PQ  determined  from  equa¬ 
tion  VI-25.  Equation  VI-46  is  identical  with  the  definition  of  a  in 
equation  III -25  if  the  tangent  moduli  and  wave  velocities  are  taken  at 
one  stress  level.  It  is  used  here  as  a  measure  of  the  hysteretic  energy 
loss. 

The  input  variables  consist  of  the  following: 

Variable  _ Description _ 

ESC  Number  of  depths  of  interest 

ZA.(N)  Specific  depths  of  interest  (ft) 

B  '  1  - 

^  +  ,  01  is  given  by  the  equation 

VI-46  (this  B  should  not 
be  confused  with  the  B  used 
in  the  body  of  the  paper) 

P  Stress- strain  curvature  param¬ 

eter  n  (see  fig.  VI-1)  •• 

Q,  Dummy  number  =  1 

EPS  Acceptable  error  level  «  0.0005 

DT  Time  increment  (sec),  a 

constant 

T21  Initiation  time  for  the  first 

unloading  wave  (sec),  t..  in 
fig.  VI-5 

TEED  Terminating  time  (sec) 

The  input  sequence  and  formats  are  as  follows: 
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Format 


Input  Card 


(A) 

[NX,  (ZA(N)  ,N=1,NX)] 

[I10,7F10.3,(8f10.3)] 

(B) 

[B,P,Q,EPS,DT,T21,TEND] 

(8fio.3) 

The  output  variables  are  given  in  three  groups  as  follows: 

Group  1.  Input  data  from  card  form  (B)  of  the  input  variables  fol¬ 
lowed  by  some  intermediate  calculations  for  debugging 
purposes. 

Group  2.  Input  data  from  card  form  (A)  of  the  input  variables. 

Group  3*  Nondimen sional  stress  (S),  time  (T) ,  particle  velocity  (V)} 
and  depth  (z),  given  by  equations  VI-42,  VI-43,  VI -44, 
and  VI-45 ,  respectively. 

A  complete  FORTRAN  II  listing  of  the  program  and  an  example  of  typical 
input  and  the  resulting  output  is  presented  in  Appendix  VI. 
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APPEEDIX  VI 

NONLUffilAR  HYSTERETIC  MODEL 
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*F  I Z 


OCT  9  /OCT  9  ,  CARD  *AAU*  FLOATING  PoINT 


C  PROGRAM  41-S8-Z0-010  29  MARCH  1968 

C 

C  NONLINFAR  HAVE  PROPAGATION  tHYSTERETIC  MODEL  1 
C 

C  CONVERTED  TO  FORTRAN  FROM  ALGOL 

C 

C 

c  *****  INPUT  paramfters  ***** 

c 

c  **  NX  NUMBER  OF  DEPTHS  OF  INTEREST 
C  **  ZA  LIST  OF  DEPTHS  OF  INTEREST  {NX  VALUES] 


C 

C  **  B  NEWNARKS  beta  parameter 

c  **  n  stress-strain  curvature  parameter 

c  **  M  M  *  N  FOR  NONSHOCK  INPUT  PRESSURE 


C  M  *  1.0  FOR  SHOCK  INPUT  PRESSURE 

C  **  EPS  ACCEPTABLE  ERROR  LEVEL 

C  **  DT  TIME  INTERVAL  BETWEEN  START  OF  UNLOADING  WAVES 
C  DURING  INITIATION 

C  **  T21  INITIATION  TIME  FOR  FIRST  UNLOADING  WAVE 
C  **  TEND  TERMINATION  TIME.  WHEN  AN  INITIATION  TIME  FOR  AN 
C  UNLOADING  WAVE  EXCEEDS  TEND,  THE  PROGRAM  ENDS 

C 
C 

C  *****  OUTPUT  PARAMETERS  ***** 

C 

C  **  S  NONDIMENSIONAL  STRESS 

C  **  T  NONDIMENSIONAL  TIME 

C  **  V  NONDIMENSIONAL  PARTICLE  VELOCITY 

C  **  Z  NONDIMENSIONAL  DEPTH 

C 


COMMON  S  1 11, 11 1, T  1 11, U I. V 1 11, 11 l.Z 111, 111 
COMMON  SAI4T l,TA  1401, VA 140 1.ZAI401 
COMMON  NJ.NX.JMAX 
C 
C 
C 

1  FORMAT! I10.7F10.3, I8F10.31! 

2  format t sfio . 3 j 

3  FORMAT {1H1.30X.22H*****  INPUT  DATA  *****,//, 9X» 1HB. 11X , 1HP, llx, 
11HO»10X»3HEPS,9X.2HDT,10X.3HT21,9X,4HTENDJ 

4  FORMATtlH0,10Fl2.5| 

1000  FORMATilHO,lP8F14.6) 

1001  FORMATI1H0.10X.18HDEPTHS  OF  INTEREST! 

C 

C 

c***  BEGIN 
C 

100  CONTINUE 
C 

C  **  READ  DEPTHS  OF  INTEREST 
C 

READ  l,NX»fZACN!,N*l,NX! 

C 

C***  AGAIN 


0001  0 
00020 
00030 
00040 
00050 
00060 
00070 
00080 
00090 
00100 
00110 
00120 
00130 
00140 
00150 
00160 
00170 
00180 
00190 
00200 
00210 
00220 
00230 
00240 
00250 
00260 
00270 
00280 
0029 
0030 
0031 

0033 

0034 

0035 

0037 

0039 

0040 

0041 

004? 

0043 

0044 


0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0054 

0055 
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c 


0056 


U1-D&-103  J8l.ll 
DO  102  Ksl.ll 
S1J..KJ  .«  0.0 
TIJ.K]  ■  O.n 


-  V-t J*-KJ — «  0.Q 
102  ZIJ.KJ  a  0.0 


C  **..  READ  AND. ,£fil NT 
C 


INPUT  PARAMETERS 


R£AH  E»A»JL»-Qx£ESjDTjTI2j1  1  .TEND 
PRINT  3 

- gfi-LNT  4.BiP.QiFPS.DT,Tt2,11.TSND 

C 

C  .**  ..COMRUTF  CONST  AMTS 
C 


-  VN  «  ■0^5-.*..A«/f g.»PI 
CN  »  0.5  -  l./t2.*P) 

- SN  a  SORTFfPI  . 

SM  «  SdRTFIQl  . ' 

BSN  e  B./.SN 

BSR  *  BSN/VN 

RM  *  2^*SM/IQtl.  I 

ALP  =  fl-B*SM/SNI/Il*B*SM/SNI 

- £81JlT__lM0^\mA£J!USN#SM.BSN1BSR.RM,ALP 

Sfl.ll  «  1.0 
..  J  =  2  .... 

Tf J.2]  a  TIJ.ll  *  2.*DT 

104  IEt  ABSFmj.21-TlJ.il 1 -DT 1 110. 110. 99 
C 

-£-*• — CQUEliXE — INTERSECTIONS  ALONG  AN  UNLOADING  HAVE 
C 

99  DO  109...K»1» J 
C 

IFtK-iaiJ05.105.106 

105  TIJ.K1  a  TI2.1J  *  CJ-2)*nT 
- 5 LIjliUL-E.JEXPF  t-.ALP*T  t  J.  1 1 / i  1  .♦ALP  J I 

Sij.ll  a  EXPFI-TlJ.in 

VIJ.K)  a  (RM-BSR1*SI1, J1**VN  +  BSR*$ l J* 1 J  **VN 
ZIJ.11  a  0. 

GO  10  109 
C 

.  106  J.F1K-J  1108. 107. 107 

107  SlJ.Kl  «  IVtJ,K-i)+BSR*SfJ»K-l]**VN)/tRSR*RMl 
VIJ/KJ  a  S  I  J. K 1 *RM 

SlJ.Kl  a  EXPFILOGFISt J.K1 1/VN! 

CU  a  tl_./BSNl*t[SlJ,K-ll*S[J»KIJ/2.1**CN 
CL  a  SM*t  ISC  J-1  .K-lUSt  J.KJl/2. 1**CN 

- ILlt-KI . BlZUrlj  K.-1 1  —  Z I  J.  K-i  1  *T  ( J.  K-l  i  *CU-T  I  J-l.  K-l  1  *CL  1  / 1  CU-CL  1 

ZIJ.K]  a  ZtJ-l.K-lUtTCJ.KI-TIJ-l.K-ll  1*CL 
GO  TO  109 
C 

108  SlJ.Kl  *  0.5*CSIJ.K-l]**vN*  SC J-1.K)**VNU  C0.5/BSRUCVI J.K-lJ- 
IV I J-l, K 1  1 

- V.U..K1  ._a_V.tJjt.K-l  1  -  BSR*ISIJ,K]-S1J,K-11**VNI 

SlJ.Kl  a  EXPFCLOGFISCJ.K1 l/VNl 

CU  a  ISN/BJ*I ISC J.K-11  ♦  SI J.KJI/2, 1**CN 

CL  a  [ SN/B 1  *  C 1st J-l. KJ  ♦  S ( J.K 1 1 /2 . 1 **CN 

TlJ.K]  a  IZI  J-l.KJ-Zl  J.K-1UTI  J.K-1UCU+TI  J-l,Kl*CLJ/(CU+CLl 


0059 

0060 

0061 

006? 

_ 0063 

0064 

0065 

0066 


0070 

0071 

007? 


0077 


0080 

0081 

008? 

0083 

0085 

0.086. 

0087 

0088 

0089 

0090 


0095 

0096 

0097 

.0098 


0104. 
01  05 
0106 
0107 
0108 


0114 
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Z  f  J, K  J  a  IT!  J.Kl-Tt  J.K-lUaCU  «■  ZtJ»K-l) 

109  PRINT  1000, J.K.SIJ.K J .Tfi. K ) ,  ¥  t  J  »  K 1  .JX J.Kl .CU.CL 


_J  _a  J+l  _  _ 

T  t  J.2 J  a  TfJ-1.2) 

_ T  TU.-1,1J . 

GO  TO  104 

_C _ 

110  JMAX  a  J-l 

- _ -JMAX _  _ _ _ 

C 

.JJICAIIL.S.UREAGE-  .INTERSECT f  ON 
C 

_ DO.  Ill  Kal.JMtx _ _ 

Sfl.K]  a  S { J, K ] 

_ ILL#JKJ_«__T!  JaK  )  . 

V  1 1#  K  J  a  V  t  J*  K ] 

111 — Z.L1jKJ  a  ZIJ/K ] 

J  a  2 


0115 

” 0T1V 

0115 

0119 

0120 

0121 

JL122 


0123 

0124. 

0125 

.0126 

0127 

Jl-129- 


C***  START 

c _ _  .... 

112  IFtTU-l.H-TENDI114.113.ll3 

.113— J...*...  J-l. . . . 

GO  TO  125 

J! _ _ 


01.29 

0.130 

0131 

0132 

0133 

JUJuL 


C  **  COMPUTE  INTERSECTIONS  ALONG  AN  UNLOADING  HAVE 

e _ _ 

114  DO  121  K«1.JMAX 

C  _ _ _ 

IFfK-11115.115,118 

_ 1-15  ST « ,  2  ..AT  1 J  - 1 .  K  *11-1 1  Jr  1..1  J  .. 

TIJ.K]  a  ST  ♦  2*EPS 

_ LEI A  B  S£  ( JJ-Jj  K  J  -ST  I -EPS  1 11  7,116.116 

116  TIJ.K)  a  ST 

_ S1.J,..K1 ...  a  EXPF.f-TtJ-.Kn 

C  a  tSN/Bt*t IStJ-l,K*l)  *  SfJ. 1)1/2. 1**CN 
- ST  aITfJ,K|lT_LJ?JLtK.f-lI*Z.IJ-l«Kfl)/CI/2.  _ _ 

117  VfJ.K]  a  VI J-l,K*ll*BSR*fSf J»K]**VN-St J-1.K*11**VN1 

..  -  ZfJ.KJ  a  0,  _ 

GO  TO  121 

C  .  . . 

118  IFfK-JMAXtl20.119.H9 

- 119-  SI  J.Kl  »_mjjLK-H.t.BSR*Sf  J.K-ll«*VNI/fBSR*RMl 

VfJ.K)  a  S  t  J. K ] *RM 

.  - SLJ.K]  »-.EXPF  ( LOGF  (SIJ.K)  J/VNJ 

CU  ■  tl./BSN!*f  tStJ.K-mStJ.K)1/2.)**CN 

. .  Gl _ a  SM*.f  ISt  J-l.KI+St  J.K1I/2.  1**CN 

TIJ.K)  a  f-Zt J,K-ll*Zt J-1 ,K)+CU*T[ J, K-l J -CL*T f J-l.kl) / tCU-CL I 

_ ZULt-F-L  2L1  J-ll.t  ,K  l.tfiL*  I T I  .i,  K I  *T  t  J-l ,  K II 

GO  TO  121 
-C. . . 

120  StJ.K)  a  0.5*ISt J.K-1)**VN*SI J-1.K*1  !**VN1  «■  [ 0 .5/BSR 1  * C V I J.K-U- 
..  _1V  LJ-l,K*lt1 

VfJ.K]  a  VIJ.K-l)-BSR*fSfJ,K)-S( J.K-1)**VN) 

SfJ.K]  a  EXPFILOGFfS t J , K 1] / VN ] 


CU  a  fSN/B)*l fSfJ.K-ll+Sf J.Kl 1/2. I **CN 
CL_«  JSN/Bj*IfSl J-l, K+ll+SI J.Kl J/2. J**CN 

Tt J.Kl  a  IZ! J-l.K+H-Z! J.K-l)*CU*Tt J,K-11*CL*Ti J-l.K+1) )/ICU*CLI 
Z  IJj  K  l_a  ZIJ,K-lJ  +  m  J.Kl-TI  J.K-1I)*CU 


0135 

013.6 

0137 

0138 

0139 

-MAIL 

0141 

0142 

0143 

0144 

-M4&. 

0147 

0149 


_ 0152. 

0154 

0155 

0156 

0157 


0163 

016A 

0165 

0166 

0168 


0173 
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c  ot76 

121  PRINT  1000»J»K#SIJ»K1,T[)»K)»VIJ»K),7(J»K)»CU,CL 
IFIJ-1U122.123.123 

122  J  =  J+l  0178 

GO  TO  112  0179 

123  NJ  «  J  0180 

C  0181 

C***  PROCEDURE  INTERPOLATE  AND 'print  0182 

C  0183 

PRINT  1001 

PRINT  1000, IZAINJ .Nsl.NXt 

CALL  INTERP  0184 

C  0188 

C  **  RESET  INDEX  TO  ZERO  0186 

C  0187 

DO  124  K=1,JMAX  0188 

Sll.KJ  ■  SIJ.K1  0189 

TI1.K1  ■  TIJ.Kl  0190 

V [ t.K  J  ■  VIJ.K)  0191 

124  ZI1.K1  =  ZlJ.Kl  0192 

J  =  2  0193 

C  0194 

C***  GO  TO  START  0198 

C  0196 

GO  TO  112  0197 

125  NJ  a  J  0198 

PRINT  1000.NX.NJ. JHAX 

CALL  INTERP  0199 

C  0200 

C**«  GO  TO  AGAIN  [NEXT  CASE1  0201 

C  0202 

GO  TO  101  0203 

CALL  EXIT  0204 

C  0205 

C***  HALT  0206 

C  0207 

END  0206 

04172.  15352 


, kp  noooooo 
S  001740*  O0«00n3 
T  0017022  00*0006 
V  0016**0  0040011 
Z  0016056  0040014 
SA  0015736  O020050 
TA  0015616  O020050 
V*  0015476  n020050 
Z*  0015356  0020050 
NJ  0015355 
NX  0015354 
JMAX  0015353 
700001  00000*0 
700002  0000050 
7U0003  0000053 
70000*  0000115 
701000  0000122 
701001  0000127 
700100  00001*2 
N  0000017 
*00001  0000021 
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In  o  I  I  |  o  I  1  o 'o  o  o 


SUBROUTINE  1NTERP 

"procedure  interpolate 


COMMON  SItl.lll.TIll.lll.Vill.nl  ,  ZJ 11 ...  11 1 _ 

COMMON  SAI40) .TAI40) .VAI40! .ZAI40) 

COMMON  NJaNX . JMAX . .  ..  . 

1  rORMATrlH0.//,21X,40H*<>**a  INTERPOLATION  ALONG  THE  HAVE  ****♦» 
1//.11HWAVE  NUMBER. 17X, IMS* 14X, 1HT, 14X, 1HV. 14X. 1HZ . / 1 
-2  FQRMATflHO.// , 21X. 39H**»**  INTERPOLATION  BETWEEN  WAVES  *«»♦». 

1//.11HUAVE  NUMBER. 17X, IMS. 14X.1HT.14X.1HV.14X.1HZ./ I 
3  FORMAT f 1H  .I6.14X.1P4E15.S./, 1 21X. 1P4E15 . 5 1  I 


PRtNT  1 
_ DO  20 B 


DO  205  N«1 » NX 
_  IF  f  N-l 1201.20 0,2  0L 

200  IF  IZA (Ml -Z[J. JMAX JJ202. 202, 208 

201  IFfZAtMl-Zt J. JMAX1 1 202,202,206 

202  DO  203  M*1 , JMAX 


0001 

00020 

00030 

.00040 

00050 

TooTo" 

_0008.fi 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

.0016.0 


IF (Z I J.M1-ZA (Nil  203, 203.204 

_ 20.3  JC.0JN.1I.NUE. _ : _ 

204  M  a  KT 

_ K_a_KI-l . . . 

ST  «  IZAtNl-ZIJ.K)  1/IZU.Ml-Zt J.KI1 

_ SA10L1  «  S  t  J.K 1  tST*,LS [ J.M.l rSI J..K11 _ 

TAtNI  a  Tt J.KUSTMTt J.M1-TIJ.K1I 
_ VAtNl  a.  VI J,  K.L*ST*  I V  [  J  ,_H  1  -  V  t  J .  K 1  I 

205  CONTINUE 
_____ ... 


00260 


206  M  a  N-l 

_2_0.2_P«U.NT_i#_Jj  IS  AJ  N  L,  TA.LNJ,  VL4.I  N  L..Z.A1  N 1 ,  N»1 ,  M  J 
20B  CONTINUE 


PRINT  2 
N.I1  a 


DO  224  Jal.NJl 

~  IF  IN-1 1 211. 209, 211 
_ 2JL9  IF  IZA  IN  1  »Z  t J, JMAX!  ]  213,213,  210 

210  IFIZAtNl-Zl J*l» JMAX I  I  213. 213.224 

211  IFfZAtNl-Zf J, JMAX] 12 13.2 1 3j,_2 12 _ 

212  IF  I ZA In] -Z l J*l» JMAX 1)213. 213.222 

213  DO  214  K«l. JMAX 
KT  ■  K 

_ IF  IZIJ.K1-ZA INI  1214,214,215 

214  CONTINUE 

2 


IFIK-JMAXt217.216.217 
_  217" IF IZ I J+l.K-1 1 -ZA IN  1121 8. 218,219 


218  KT  ■  J*1 

H  a  K-1 _ 

ST  ■  t7AfN)-ZtKT#MJ 1 / 1 Z l j#  K 1 “Z  t  KT#  M  J | 
SAIN]  a  SlKT#M)*ST*t$l J#Kl-StKT,M|) 

TAfN!  a  TlKT#MJ*ST*ITl J.K) -T tKT.M) ) 

- StAJLNJ _ VIKT#M1*ST*(VI  J.KJ-VIKT.Ml  1 

GO  TO  221 

_ 219  K  »  K-i _ 

220  M  =  K 

...S.T-»  JZAtNl-Zt J.KH/tZI J«1»H1-Zt J.KI] 

SA  t  N ]  a  St J»K1*ST* tS{ J+1.M1-S1 J,K1 1 

— . -TAJNL.a.  II J#Kl*ST*tTl J+l.Ml-Tt J#K)  J 

VAtNt  a  Vt J.K)*ST*tVt J+l.MJ-Vt J.KJ1 

_ 221_.CQ.NT  I  Nli£  . 

C 

. . N  _a  .  N_ 

GO  TO  ?23 

222  . M  a  N-l_  .  . 

223  PRINT  3.J*tSAtNJ»TA[N),VAtN)»ZAtN1#N«l.M1 

—  224  COMT  t  Ml  IF _ 

C 

.  _ RETURN .  .. 

END 

_  .O.iAZfL#..  .15.352  - . 


INTERP  0000012 _ 

S  0 01T4q4  p0^oon3 
J .  0017022  Q 040014 
V  0016440  0040017 

Z _ 0 016056  00 4 0 0_2 2 

SA  0015736  0020050 
-I* — nniSftife  nn2nn5n 
VA  0015476  0020050 


NJ  0015355 
NX  nni5354 
JMAx  0015353 
■/■Doom  nnnno45 
/00002  0000112 

/DO 00.3  00  0  015.6 _ 

J  0000025 


N  0000027 

/tio2ni  nnnt)23n 

'00200  0000210 


'002o8  0000572 

-'0020-6— omn535 _ 

M  0000032 

-SI — nano  o  3.3 _ 

'002p3  0000275 

-'002q4  ..00011104 _ 

K  0000034 

JSI - go  00036 _ 

'002o5  0000523 

'002n7  nnnn54n _ 

NJ1  0000350 

-'.00211  Q.QQQ662 _ 

'00209  0000621 
-'.0021 3  nnnn723 _ 


0067 


0074 

0075 

0076 
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/00101  0000165 
J  0000.022 
-00032  0000023 
K  0000024 
0.0  0000026 
'00102  0000224 
'00103  0000243 
B  00000.3.0 
P  0000032 
0  0000034 

BPS  0000?62 
DT  0000264 
♦C0002  0000266 

TEND  . O00fi2?o. 

VN  0000272 
0-5  0000274 
1*  0000276 

2*  0000300 

CN  0000304 
SN  .  .0000310 

SQRTF  0140410  EXT  PROG 

SM  0000412 

BSN  0000414 

BSR  0000416 

RM  0000420 

ALP  0000422 

1*0  0000430 

/00104  0000576 

ABSF  0140432  ExT  PROG 

/00110  0002206 

/0 0 099  0000623 

/001C5  0000633 

/00106  0001045 

EXPF  0140434  EXT  PROG 

. I R  C140710  EXT  PROG 

0.  0000712 

/00109  0002061 

'OOlO8  0 0fll4l4 

/00107  0001053 

LOGF  0140714  EXT  PROG 

CU  0000716 

CL  0000720 

/00111  0002271 

'00112  0002321 

'00114  0002343 

'00113  0002337 

'00125  on04122 

'00115  0002353 

'00118  ono2655 

ST  000.1516 

'00117  0002565 

'00116  0002435 

C  0001520 

'00121  0003666 

/00120  0003250 

'00119  -0J102663. 

♦00011  0003164 
'00122  0003756 
'00123  0003762 
INTERP  0103166  EXT  PROG 
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/00124  0n0407l 

EXIT  0104140  EXT  PROG 
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0.  04 

0.06 

0.03  0.1 

0.2 
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ZA  (2) 
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8 


input  data 

EPS 


p 

0.64667  2.00000 

7.5DC080E-01  2.5*000*E-01 

2.o*o*oc6i6o '  1. *008**6*00 
2.0001006.00  2.*0000QE*00 

3.000-006*00  1 • 00000*6*00 

3.000*0*6*00  2-00000*6*00 

3.00**00E*00  3. *OOO0OE*O* 
4.0001006*00  1.000000E»00 

4.ooo*ooeTob  "2.~*obooo6*oo 
4.000*006*00  3.0000006*00 

4.000*006*00  4, *0000*6*00 

5.6oo*OOE*O0  1 <0000002*00 
5.000*00E*00  2.00000*6*00 
5. 000*006*0*  3.1000006*00 

5.000*006*00  A .  0OOOOOE«OO 
5.000100e*00  5.ooooooe*oo 
2-OOOlOOE+oo  1  *0000006*00 
2.000*006*00  2.10000*6*00 
2.ooo*oce*oo  3.1000006*00 
2.000*006*0*  A.  t(l0O00E*0O 
2.i)no*00E*0*  5.0000006*00 

3.000*0*6*00  1 . 0O000OE*00 

3, 000*006*00  2, 0000006*00 

3.000*006*00  3. *000006*00 

3.000*006*00  A. *000006*00 
3.000100E*00  5. *000002*00 

4.000*006*00  1 .0000006*00 

4.000*006*00  2*0000006*00 

4.000*006*00  3.0000006*00 


S 

1.00000 
1.4142146*00 
0.0000056-01 
0.0016016-01 
0.0060106-01 
0.0073736-01 
0.0083566-01 
0.0820166-01 
0.0030676-01 
0.0840506-01 
0.0850326-01 
0.0780246-01 
0. 

0.0797456-01 
0.0807276-01 
0.0817006-01 
0.0721076-01 
0.0731706-01 
0.0741616-01 
0. 07514J6-01 
0.0766406-01 
0.0610046-01  3.808102E-03 

0.0620766-01  4.4553116*03 

0.0630576-01  5.0124126*03 

0.0645556-01  5.826120E-03 

0.0673116-01  7.3778406*03 

0.0500016-01  5.0124306*03 

0.0500826-01  5.569647E-03 

0.0524706-01  6.3834776-03 


TEND 
15.00000 
6.285304E-01 
0. 0000006-01 
1.801820E-03 
0. 0000006-01 
1.4675846*03 
2.6485546-03 
0.0000006*01 
1.0434156-03 
2.2243166-03 
3.405701E-03 
O.OOOoOOE-Ol 
6.1031226-04 
1.8001456-03 
2.0810516*03 
4.162n326-03 
0.0000006*01 
1.1807506*03 
2.3615526-03 
3.5424286-03 
5.2*73026-03 
0.0000006-01 
0.0770446*01  1.1*0*016*03 

0.0774816*01  2.3613736*03 

9.976774E.01  4.0850636*03 

9.075473E.01  7.3748326*03 

0 • 072246E-01  0. 0000006-01 

O.O7J7B2E-01  1.1805716*03 

O.O71076E-01  2.0048786-03 


1 .0000006*00  3.50245*6-01 

-6.2*07666-40  -6.28076*6-40 
2.1208356*00  O.OO8O60E-OI 
2.1208356*00  9.99896bE*01 

2.1206146*00  2.1207686*00 
2.1206766*00  0.0075046-01 
2.1206766*00  0.0075046-01 
2.1203046*00  2.12053*E*00 
2.1204486*00  2.120588E.00 
2.1205006*00  0.0066726-01 

2-1205006*0(1  0.0066776*01 

2.1201746*00  2.120307E*00 
2.1707106*00  2.12036*6*00 

2.1702716*00  2.1204176*00 

2.1703246*00  0.00584*6-01 

2-1203246*00  0.00584*6*01 

2.11 O870E*00  2.12007iE»00 

2.1109236*00  2.120123E-00 
2.1100756*00  2.1201786*00 
2-1200416*00  0.00470*6*01 

2.1200416*00  0.00470*6*01 

2.1102806*00  2.11962*6*00 

2.1103326*00  2.110680E*00 

2.1193086*00  2.1197506*00 

2.1105116*00  0.0920876*01 

2.1105116*00  0.9920876-01 

2.1186006*00  2.1lO037E*00 

2.1187566*00  2.1101176*00 


»T  T21 

0.00010  0.00040  0.00100 

I.O1OOOOE*O0  4.7140456-01 
1.0*00006*03  9.0045516*01 

1.8020176*03  9.0037616*01 

1.4"0C00E*0S  0.0923736*01 

2.0020566-03  O.OO1730E-OI 
2.6480406-03  9.001266E-01 

1.B10000E-03  0.0001056*01 
2.2020856-03  0.080600E-01 

2.64BO06E-0S  9.089236E-01 

3.4*50206-03  0.088772E-01 

2.2*00006-03  O.O880I8E-01 
0787636-01  2.4021046-03  0.0876606*01 

3.0400436*03  0.0872066*01 
3.6*60036*03  0.0867426*01 

4.1620856*03  0.0862706*01 

2.7842006*03  0.0845726*01 

3.3412016*03  0.0841006*01 
3.8082036-03  0.0836456*01 

4.4552276*03  0.083182E-01 
5.2688316*03  9.0824756*01 
0.0784076*01 
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4.00onOOE*00 

4.aoOOOOE*OQ 

9.95523SE-01 

7.935404E-03 

-9.969775E-01 

8.1937076-03 

2.110869E+OO 

2.11*2646*00 

4.00nnOOE«00 

5.noooooE*on 

9.957989E-01 

9.4878546-03 

9.96B475E-01 

9.482*736-03 

2*11 9oi6E*0  0 

9.9906496-01 

5.0000006*00 

i.nofloooE*oft 

9.938919E-01 

6.1260656-03 

9.9660**6-01 

0.0(100006-01 

2.1190166*00 

9.9906406-01 

5.  oooto06«oo 

2.ooo6ooe*oo 

9.94n415£-01 

6.940847E-03 

9.965381E-01 

1.724105E-03 

2.1181136*00 

2.1184786*00 

5.OOO''O0E*OO 

3.nonooos«-op 

9.943170E-01 

8.4920896*03 

9.964081E-01 

5.011895E-03 

2. 1182276*00 

2.1186216*00 

S.OOP'XIOE+OO 

4,0000006*00 

9.945924E-01 

1.004566E-02 

9.962701E.O1 

8.3flln2lE«03 

2 . 116373F«  00 

2.11876*6*00 

5/oonoooE*oo 

5.100000E*00 

9.948675E-0t 

1.1598776-02 

9.9614826-01 

‘iri59i32E-02 

2.1185206*00 

9.988313E-01 

6.0000006*00 

1.0000006*00 

9.9227526-01 

7.754828E.03 

9.957040E-01 

0.0000006-01 

2.118520E*00 

9.9883136-01 

6. 000°COE*00 

2.0(100006*00 

9.925506E-01 

9.3073706-03 

9.955740E-Q1 

3.2871756-05 

2-1172856.00 

2-1170302*00 

6.000000^00 

J-6oooooE*oo 

9.928258E-01 

1.0060356-02 

9.994440E-01 

6. 37591 Of* 03 

2-1174326*00 

2-1179766*00 

6.000^00E400 

4.0000006*00 

9.931000E-01 

1.241376E-02 

9.9531416-01 

9.88501*6-03 

2.1175796*00 

2.1101236*00 

6»000f»D0E*00 

*-inooooE‘oo 

9.935030E-oi 

1.4683346-02 

9.951239e-C1 

1.4f>7^4gE*02 

2‘.H7759E*od 

9.9054376-01 

7.00nftOOE«00 

1  aO(IPOOOE*00 

9.8919B8E-01 

1.0850916-02 

9.93990*6-01 

o-ooonobE-oi 

2.11775*6*00 

9.9854396-01 

7.{l00'»00E+00 

2.0(100006.00 

9.894739E-01 

1.241380E-02 

9.93860*6-01 

3.2874686-03 

2.1156426*00 

2.1l66l?E*00 

7.0001006*00 

3.000000E*00 

9.8974B7E-01 

1.3967836*02 

9.937305e-01 

6.5754716-03 

2.1157896*00 

2.1167586*00 

7.0001006*00 

4.'»00003E*OP 

9.9015136-01 

1.6238276-02 

9.935402E-01 

1.137966E-0? 
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VII 


Comparative  Studies  of  Various 
Stress-Strain  Models 

A.  Purpose 

A  one -dimensional  -wave  propagation  problem  in  a  typical  soil  is 
solved  herein  utilizing  the  computer  programs  presented  in  Chapter  III, 

V  and  VI  for  the  linear  hysteretic  model,  the  nonlinear  locking  model  and 
the  nonlinear  hysteretic  model  respectively,  in  order  to  compare  the  re¬ 
sulting  solutions  with  a  finite  difference  solution  (28)  of  the  same 
problem  in  which  the  "actual"  stress-strain  curve  for  the  soil  is  used. 

The  stress-strain  relation  for  the  soil  under  consideration  is  shown 
in  fig.  VII-1.  During  the  early  stages  of  loading  the  stress-strain 
curve  is  concave  downward  indicating  a  "plastic"  behavior  of  the  material. 
However,  as  the  stress  increases  further,  the  curve  becomes  concave  up¬ 
ward  indicating  a  "strain -hardening"  behavior  of  the  material.  At  stress 
levels  higher  than  those  indicated  in  fig.  VII-1,  the  tangent  modulus  will 
generally  continue  to  increase.  During  unloading  the  soil  exhibits  an 
extremely  stiff  behavior  at  all  stress  levels  except  near  zero  where  a 
sharp  breaking  tail  is  evident.  The  result  of  a  complete  load  cycle  is  a 
considerably  large  permanent  strain. 

The  overpressure -time  relation  used  for  the  comparative  studies  is 
shown  in  fig.  VII-2  and  is  described  by 

P(t)  =  700  e-t/38  VII-1 

where  P(t)  is  expressed  in  psi  and  t  is  expressed  in  msec.  The  initial 

O 

mass  density  of  the  soil  is  P  =  3.66  slugs/ft'5  . 
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Overpressure  (psi) 


Fig.  VII-2.  Overpressure-time  relation  for  model 
comparison  study. 
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B.  Approximations  to  the  soil 

stress -strain  curve  with  the  various  models 

The  approximations  to  the  stress-strain  curve  of  fig.  VII-1  for  the 
various  models  are  shown  in  figs.  VII-3  through  VII-6.  Fig.  VII-3  depicts 
the  piecewise-linear  fit  which  was  used  to  obtain  the  finite  difference 
solution  of  the  problem.  Six  straight  line  segments  were  used  to  represent 
the  loading  stress-strain  characteristics  and  four  straight  line  segments 
were  used  to  represent  the  unloading/reloading  characteristics.  The  com¬ 
parison  (shown  in  fig.  VII-3)  of  the  piecewise-linear  fit  with  the  actual 
soil  relation  indicates  that  the  piecewise-linear  model  does  indeed  accu¬ 
rately  represent  the  "actual"  soil  stress-strain  relation. 

Fig.  VII -4  shows  the  nonlinear  locking  model  approximation  to  the 
soil  stress-strain  curve  (this  model  was  discussed  in  Chapter  V)  and  figs. 
VII-5  and  VII-6  show  the  linear  hysteretic  model  and  the  nonlinear  hyster- 
etic  model  approximations,  respectively  (these  models  were  discussed  in 
Chapters  III  and  VI).  For  ease  of  comparison,  the  actual  soil  stress- 
strain  relation  is  also  shown  (as  a  dashed  line)  on  each  of  these  figs. 

Note  that  all  of  the  approximations  to  the  loading  stress-strain  relation 
pass  in  common  through  two  points  on  the  actual  curve,  the  zero  stress 
level  point  and  the  700-psi  stress  level  point. 

C .  Results  of  the  computer  runs 

The  results  of  the  computer  runs  for  the  various  models  are  presented 
as  plots  of  the  attenuation  of  peak  vertical  stress  and  particle  velocity 
with  depth  in  figs.  VII-7  and  VII-8.  In  addition,  stress-time  histories 
at  various  depths  are  shown  in  fig.  VII -9  for  the  linear  hysteretic  model 
and  in  figs.  VII -10  and  VTI -11  for  the  piecewise-linear  model.  Particle 


Vertical  Stress  (psi ) 


Fig.  VII-3.  Piecewise-linear  stress-strain  fit  for 
use  with  finite  difference  solution. 
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Vertical  Stress  (psi) 


Vertical  Stress  (psi ) 


Fig.  VII- 5 .  Linear  hysteretic  model  fit. 
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Vertical  Stress  (psi) 
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velocity-time  histories  at  various  depths  are  shown  in  figs.  VII-12  and 
VII-13  for  the  piecewise-linear  model. 

D .  Discussion  of  results 

The  attenuations  of  peak  vertical  stress  and  particle  velocity 
■with  depth  (figs.  VII -7  and  VII-8,  respectively)  for  the  various  models 
clearly  indicate  that  the  attenuation  rates  are  related  to  the  hysteretic 
energy  loss  potential  (area  under  the  stress-strain  curve)  as  well  as  the 
details  of  the  stress-strain  curve  associated  with  each  model.  The 
piecewise -linear  model  (fig.  VII -3)  provides  the  highest  attenuation  rate 
of  peak  vertical  stress  and  particle  velocity  with  depth.  This  is  partly 
due  to  the  fact  that  the  plastic  behavior  included  in  this  model  at  the 
lower  stress  levels  causes  a  spreading  and  a  slowing  of  the  wave  front 
(i.e.,  an  increase  in  rise -time  and  a  decrease  in  wave  velocity)  as  the 
wave  propagates  down  the  soil  column  (see  figs.  VTI-10  through  VII-13), 
and  partly  because  the  piecewise-linear  model  incorporates  a  high  hys¬ 
teretic  energy  loss  potential.  The  spreading  and  slowing  of  the  front  is 
unique  to  this  model;  in  effect,  these  phenomena  allow  unloading  waves  to 
have  more  influence,  i.e.  cause  more  attenuation,  at  a  given  depth. 

The  other  models  all  produce  shocks.  A  shock  will  develop  if  the  entire 
loading  stress-strain  relation  is  of  the  form 


d2g 

de£ 


s  o 


(which  is  the  case  for  the  nonlinear  locking  and  the  linear  and  nonlinear 
hysteretic  models). 

The  nonlinear  hysteretic  model  (fig.  VII-6),  having  the  least  amount 
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Depth  (ft ) 


Peak  Vertical  Stress  (psi) 
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Peak  Vertical  Particle  Velocity  (fps) 

0  10  20  30  40 


Fig.  VII-8.  Attenuation  of  peak  vertical 
particle  velocity  with  depth;  various 

models . 
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of  hysteretic  energy  loss  potential,  results  in  the  lowest  attenuation 
rates.  The  linear  hysteretic  (fig.  VII -5)  and  the  nonlinear  locking 
(fig.  VII -4)  models,  having  practically  the  same  hysteretic  energy  loss 
potentials,  result  in  similar  attenuation  rates  of  the  peak  stress  and 
particle  velocity. 

E.  Other  ground  motion  parameters 

An  interesting  parameter  associated  with  wave  propagation  (not  ex¬ 
plicitly  studied  herein),  as  far  as  the  details  of  the  stress-strain  curve 
are  concerned,  is  particle  acceleration.  For  a  shock-fronted  overpressure 
pulse  (fig.  VII -2)  the  stress  and  particle  velocity  wave  forms  will  retain 
a  shock  front  if  the  loading  stress-strain  properties  of  the  material  are 
of  the  form  d  o/de  2  o  ,  as  was  discussed  in  Chapter  II.  This  "jump" 
discontinuity  at  the  wave  front  results  in  an  infinite  value  of  particle 
acceleration .  Therefore ,  for  the  nonlinear  locking ,  nonlinear  hysteretic 
and  linear  hysteretic  models,  the  peak  acceleration  will  he  infinite  at 
all  depths.  However,  the  piecewise -linear  model  cannot  sustain  a  purely 
shock -fronted  particle  velocity  or  stress  wave  form  due  to  the  inclusion 
of  the  plastic  behavior  of  the  loading  stress-strain  curve  at  the  lower 
stress  levels  (see  Chapter  II).  This  is  clearly  demonstrated  in  figs.  VII- 
10  through  VII -13  where  the  rise  times  of  the  stress  and  particle  velocity 
wave  forms  can  be  seen  to  continuously  increase  with  depth.  This  results 
in  finite  acceleration  values  and  a  large  attenuation  rate  for  peak  accel¬ 
eration.  Unfortunately,  accelerations  obtained  from  finite  difference 
solutions  have  only  qualitative  value  (28);  hence  no  computed  accelerations 
are  presented  here. 
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Fig.  VII-9 .  Typical  stress-time  histories  for  several  depths;  linear 

hysteretic  model. 
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Fig.  VII-10.  Typical  stress-time  histories  for  several 
depths;  piecewise-linear  model. 


Fig.  VII-11.  Typical  stress-time  histories  for 
several  depths;  piecewise-linear  model. 
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Fig.  VII -12.  Typical  particle  velocity-time  histories  for 
several  depths;  piecewise-linear  model. 
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Fig.  VII-13.  Typical  particle  velocity- time  histories  for 
several  depths;  piecewise-linear  model. 


Particle  displacements  are  not  discussed  herein  either  because  the 
computer  codes  presented  in  Chapters  III,  V  and  VI  do  not  calculate  dis¬ 
placements. 

F .  Conclusions 

Figs.  VII -7  and  VII -8  clearly  indicate  that  computed  -wave  propagation 
results  are  sensitive  to  the  stress-strain  model  chosen  to  represent 
actual  laboratory  data.  However,  the  variations  in  results  exhibited  by 
the  models  considered  for  this  study  are  not  highly  significant;  in  fact, 
the  "scatter"  of  these  attenuation  curves  is  of  the  same  order  of  magnitude 
one  would  expect  to  see  in  data  retrieved  from  actual  laboratory  one- 
dimensional  wave  propagation  tests  (29). 
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